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1.   INTRODUCTION 

Optimization  techniques  can  be  divided  into  two  classes,  single  stage 
and  multistage.   In  multistage  optimization  techniques,  a  certain  relation- 
ship is  used  to  isolate  the  interconnections  between  the  various  stages. 
Thus  one  stage  is  searched  at  a  time  instead  of  all  the  N  stages  simul- 
taneously.  In  this  way,  an  N-dimensional  problem  is  converted  into  N  one- 
dimensional  problems  if  the  problem  has  only  one  control  variable.   The 

« 

multistage  optimization  techniques  can  be  classified  into  classical  tech- 
niques (calculus  of  variation)  and  dynamic  programming. 

In  case  of  calculus  of  variation,  the  resulting  equations  form  a  two- 
point  boundary  value  problem  (2,15).   The  differential  equations  encountered 
in  practical  applications  are  generally  nonlinear  and  cannot  be  solved  an- 
alytically.  Finding  numerical  answers  for  this  nonlinear  boundary  value 
problem  is  very  tedious  especially  if  there  is  a  large  number  of  equations 
with  a  large  number  of  initial  values  missing.   This  has  limited  the  use 
of  the  calculus  of  variation. 

The  maximum  principle  is  a  very  powerful  tool  for  obtaining  analytical 
solutions  of  linear  optimization  problems  with  inequality  constraints  on 
control  variable  (7) .   But  when  the  problem  is  nonlinear  and  an  analytical 
solution  cannot  be  obtained,  the  maximum  principle  gives  rise  to  similar 
boundary  value  difficulties. 

Dynamic  programming,  although  free  from  the  boundary  value  difficulty, 
has  a  serious  drawback  because  of  its  storage  requirements  on  the  computer. 
Instead  of  solving  any  individual  process,  the  dynamic  programming  technique 
solves  a  family  of  related  processes  (20).   In  here,  as  in  the  other  multi- 
stage techniques,  the  problem  of  an  N  dimensional  search  is  reduced  to  N 


.sional  search  problems  if  the  problem  has  only  one  control  vari- 
However,  in  Investigating  one  stage  at  a  time,  all  possible  combin- 
ations of  the  stage  variables  for  the  previously  calculated  stage  must 
be  stored  in  the  memory  of  the  computer. 

This  storage  requirement,  often  referred  to  as  the  "curse  of  di- 
mensionality," becomes  too  excessive  to  permit  the  use  of  dynamic  program- 
ming for  a  problem  in  which  more  than  three  state  variables  are  involved. 
Thus,  if  a  three-dimensional  problem,  i.e.  involving  three  state  variables, 
is  to  be  solved  and  if  it  is  decided  to  have  each  state  variable  dis- 

cretized  into  50  values,  then  because  of  the  interpolation  required  in  the 

3 
dynamic  programming  approach,  (50)   values  have  to  be  stored.   Thus  it  is 

frequently  impossible  to  handle  even  a  three-state  variable  problem  with 
straight  forward  dynamic  programming. 

Thus  it  is  seen  that  the  dimensionality  difficulty  in  dynamic  program- 
ming and  the  boundary  value  problem  in  the  classical  methods  limit  the 
number  of  state  variables  in  a  problem  that  can  be  treated  by  these  tech- 
niques.  It  should  be  noted,  however,  that  these  two  difficulties  are 
totally  different  from  each  other.   The  dimensionality  difficulty  requires 
more  computer  memory  while  the  boundary  value  demands  more  computer  time. 
Also,  the  classical  boundary  value  problem  approach  represents  an  iterative 
procedure  to  obtain  the  numerical  solution  while  dynamic  programming 
represents  an  expansion  of  the  original  problem. 


2.   GRADIENT  TECHNIQUES 

The  methods  of  gradients  seem  to  remove  the  difficulties  experienced 
in  dynamic  programming  and  the  classical  multistage  techniques.   Although 
there  are  various  approaches  with  these  methods,  the  basic  philosophy 
remains  the  same.   When  use  of  the  gradient  methods  is  contemplated,  the 
problem  is  formulated  as  a  final  value  problem.   In  other  words,  the  per- 
formance index  or  the  objective  function  is  selected  as  the  value  of  some 
function  at  the  end  of  the  process.   This  is  not  a  serious  restriction. 
Thus  if  the  performance  index  is 

Cf 
J  =  /   f(x)dt 
0 

Then 


f=f« 


Introducing  an  additional  state  variable  x  _ 


dx  +n 

n+1    _,  . 

Tt — =  f(n) 

and  Xn+l(t0}  =  °- 

The  original  integral  performance  criterion  is  replaced  by  a  criterion 
which  calls  for  extremizing  the  final  value  of  an  element  of  the  state 
vector.   Philosophically  at  least,  extremization  of  any  performance  cri- 
terion should  be  possible  by  using  the  following  approach  underlying  the 
methods  of  gradients. 

First  a  sequence  of  values  of  control  vector  is  taken.   Then  a  compu- 
tation is  made  of  the  gradient  of  the  performance  index  with  respect  to  each 


control  vector.   Next  each  control  vector  is  improved  by  moving  it  in  the 
itc  direction  alon>;  the  individual  gradients.   This  improved 
[uence  oi"  control  vectors  then  becomes  the  basis  for  the  next  iteration. 
In  the  following  sections,  the  first  variation  method,  a  technique 
suitable  for  optimizing  nonlinear  complex  problems,  is  summarized.   Then 
the  second  variation  method,  which  is  more  sophisticated  than  the  first 
variation  method,  is  discussed.   Three  applications  of  this  method  in  the 
field  of  production  planning  and  control  illustrate  the  advantages  and 
disadvantages  of  this  method. 


2.1  The  First  Variation  Method 

Because  of  its  computational  appeal,  various  versions  of  the  gradient 
methods  have  been  developed  for  optimization  calculations.   A  gradient 
technique  for  the  numerical  solution  of  dynamic  optimization  problems  is 
generally  known  as  the  functional  or  serial  gradient  technique.   This 
technique  has  been  applied  successfully  to  solve  problems  in  aerospace, 
control  and  chemical  engineering  systems  (5,6,10,16,17,20,21).   The  con- 
tinuous  version  of  the  functional  gradient  technique  was  developed  inde- 
pendently by  Kelley(lO)  and  by  Brayson  and  his  coworkers  (5) .   A  compre- 
hensive treatment  of  this  technique  and  of  the  gradient  methods  in  general 
can  be  found  in  the  article  by  Kelley(21) . 

In  this  method,  the  convergence  is  generally  independent  of  the  initial 
guess  used  in  the  iterative  procedure,  although  the  rate  of  convergence  or, 
alternatively,  the  computer  time,  is  affected  by  the  initial  guess.   The 
number  of  equations  to  be  integrated  in  the  forward  direction  is  (n+1) ; 
i.e.  these  equations  are  integrated  from  t=0  to  t=tf.   There  are  (n+1) 
recursive  equations.   There  are,  however,  no  equations  to  be  integrated 
in  the  backward  direction  from  t=t   to  t=0.   The  first  variation  equations 
are  simpler  than  those  of  the  second  variation  method. 

The  main  drawback  of  the  first  variation  method  is  that  a  very  large 
number  of  iterations  must  be  made  in  order  to  approach  the  optimal  tra- 
jectory.  More  important  is  the  fact  that  the  trajectory  approaches  the 
optimum  but  does  not  actually  reach  it  within  a  finite  number  of  iterations. 
In  some  cases,  the  trajectory  is  far  from  the  optimum  after  a  large  number 
of  iterations  and  the  rate  of  convergence  becomes  too  slow  to  permit  further 
iterations.   This  method  cannot  conveniently  handle  the  problems  with  in- 
equality constraints  on  the  state  variables. 


ond  Variation  Method 

The  pioneer  work  in  the  area  of  second  variation  method  has  been 
d  out  by  Hryson  and  his  coworkers  (4,5),  Kelley  and  his  coworkers 
(10,11),  Merriam  (25)  and  Jaswinski  (9).   Mitter  (26)  and  Breakwell  and 
Ho  (8)  have  also  added  to  the  work  in  this  field. 

This  method  is  a  natural  evolution  of  the  first  order  linearizations 
used  in  the  first  variation  method  in  which  the  equations  are  linearized 
by  truncating  after  all  linear  terms.   The  second  order  and  higher  order 
terms  are  thus  ignored.   It  is  well-known  that  the  use  of  a  linear  ap- 
proximation in  a  gradient  search  procedure  is  an  excellent  means  for  ar- 
riving near  the  optimum  point  quickly  and  from  almost  and  stationary 
starting  point.   Near  the  optimum,  however,  the  linear  approximation 
becomes  deficient  and  it  is  necessary  to  turn  to  a  second  order  approxi- 
mation to  achieve  the  optimum.   A  useful  optimization  procedure  is  to 
initially  use  the  first  variation  to  get  near  the  optimum  trajectory  and 
then  to  switch  to  the  second  order  method  for  refinement. 


2.3  Derivation  of  the  Second  Variation  Method 

Consider  a  process  which  can  be  represented  by 

^-  f[x(t),  8(t)]  (1) 

where  x  is  n  dimensional  state  vector,  0  is  r  dimensional  control  vector 
and  x(Q)    is  prescribed.   No  terminal  constraints  are  to  be  imposed  on 
_x(t, )  }  although  the  final  time,  tf,  may  be  specified.  • 

Suppose  it  is  desired  to  minimize  the  following  performance  index: 

t 

1  =  T  =  i 


f 

l[x(0),tf]  =  I  =  !     J(x,9_,t)dt  (2) 

0 


From  Equation  2,  this  equation  results 


^   »  J(x,£,t)  (2A) 


Since  the  performance  index  as  given  by  Equation  2  is  subject  to  the 
system  constraints  of  Equation  1,  consider  the  minimization  of  the  un- 
constrained performance  index  as 

*  tf       d* 

1=1  +  /   z'(f  -  -7")dt  (3) 

o dt 

where  z  is  a  vector  of  n  Lagrangin  multipliers.   Substituting  Equation  2 
into  Equation  3  results  in 

*  tf  d2L 

I   =  /   [J(x,6_,t)  +  z'(f  -  ^]dt  (4) 

In  order  to  minimize  I  ,  an  iteration  algorithm  can  be  constructed 

such  that 


tf  (j+1) 

x*(j+D  _  i       [jCj+D  +  z,(J+i)ff(j+i)   ( 

0   I        " 


dt 


dt 


(5) 


converges  in  a  desirable  way.   The  superscript  (j+1)  is  used  to  indicate 
the  number  of  iteration,  and  it  is  desired  to  have 


i*(0)  >  i*(1)  >  ...  >  i*«>  >  i*"+1)  >  ... 


(6) 


To  construct  the  desired  iterative  algorithm,  the  values  of  the 
functions  at  iteration  (j+1)  can  be  expressed  in  terms  of  the  j    iteration 
by  means  of  Taylor's  series  expansion.   Retaining  only  the  terms  up  to  the 
second  order  gives 


3xCj;  39U; 


+  i  6x(:)'lVil 

+  2  6^       (j)2  "- 

oX 


6x 


36(J).  3x(j)   " 


+  i6e(^'^ 


•  x,  a2T(j) 


2  - 


36 


(jT26^ 


(j) 


(7) 


where, 


6x<J>  =  x"+1>  -  ,0) 


66 


(j)    :   fl(j  +  1)   _   fl(J) 


(8) 


afj 

3x2 


a3i 

3x? 


32J 


82J 
axl3x2 


3x  Sx., 
n  1 


32J 

3x..3x 
1  n 


32J 


3x 


n 


a*3_ 

36  3x 


,2j 


dQ    dx 


32J 


36  3x. 
r  1 


32J 


39.,  3x 
1  n 


32J 


39  3x 
r  n 


(9) 


The  superscript  (j)  has  been  omitted  in  Equation  (9)  for  clarity, 
Thus  it  is  seen  that 


6'  J^-6x  = 


n    r    32J 


-  3  9  3x 


69.  6x. 


/.  38,3x,    i    i 
i=l  j=l   j   i 


(10) 


Next,  define  the  Hamiltonian 


H  =  z*  f 


(ID 


and  expand  H  at  the  (j+1)    iteration  up  to  the  second  order  terms  as  a 
function  of  H  at  the  j    iteration.   Note  that  H  is  a  function  of  x,  9_, 


and  _z  and  that  — r  ■  0 
3z 


10 


3z 


3x 


(12) 


(j)2  - 


e(j)3x(:i)  " 


39 


2r7(j) 


i   ,2-(j) 


+  {eo),4H^5z(J)  +  ^,4 


e(j)3z(j)  " 


3x(j).  3z(J} 


6z 


(j) 


Now  consider  the  nonlinear  performance  equations.  If  these  equations 
are  linearized  by  Taylor-series  expansions  and  by  retaining  only  the  first 
order  terms,  the,  result  is 


i 


Xi)\ 


V=— 7-  (5fCj)'/3x(j))'  6x(j)  +  (M^'/M.0')'  «i0)  d3) 


with  6x_(0)  =  0^  since  the  initial  conditions  are  constant.   This  last 
equation  may  be  rearranged  by  noting  that 


dx 
dT 


(J+l)     dx(i) 


=  6(-^-)  +  f(j) 


(14) 


11 


Thus   Equation  13   can  be  rewritten  as 


dx 


P(j+D 


dt 


(j)    +      32H(j) 


3z(j)    3x(^ 


6x(^>    + 


32H(J> 


a«<J>ae<J>    6i 


(j) 


(15) 


Furthermore, 


,(J+i)  „  z(i)  +  p(J)  6x(J) 


so  that  6z(j)  =  P(:j)  6x(j) 


(16) 


where  the  matrix  P  is  defined  by 


P  = 


3z. 

3x! 


3z 

i 

3x. 


!5 

3x 


3z 
2 

3x 


3z'  ' 


(17) 


It  is  a  symmetrical  matrix, 


i.e. 


3z  .    3z . 
i  =  _J. 

3x. 


3x 


Clearly  P_  is  unknown  explicitly  at  this  point.   For  the  sake  of 
clarity,  the  superscript  (j)  is  omitted  in  the  subsequent  derivation. 

If  now  the  normal  Hamiltonian  function  is  defined  as  H  =  J  +  z '  f_, 
then  the  above  expressions  can  be  substituted  into  Equation  4  to  yield 


12 


Cf 


1  -   I      +  /        {feJ      ^x  +    (-)      60.  +   ^   «x»  -j  fix 

0  *"        —  —  ^v 


3x 


2  2  9f ' 

+  «!•  ^  fii  +  i  «if  ^i  «I+  fie'  3^- p  fix  (18) 


91'  h  dX- 

+  fix'   t —  P   fix  -   z'fi  ^*-  -  fix'    P   6  3-  1-dt 
—     8x     —    —       —       dt  —    —       dt 


To  further  simplify  Equation  18,  use  of  the  adjoint  equation  is 

made.   Thus 


This  is  easily  obtained  by  defining 

M  =  min   [J  (x,_9,t)  +  z'f]  (20) 

9 

where  the  adjoint  variable  z_   is  defined  by 

*  =  lx"  (21) 

ox 

But  from  the  principle  of  optimality  in  dynamic  programming,  it 
follows  that  for 


Cf 
I°(x,  t)  =  min  /  J(x,9_,X)  dX  (22) 

6   0 


that 


1    (x,t)    ■   nun 


L+At  f 

/   J(>£,0_,A)dA    +      /      J(x,_0,A)dA 
t  t+At 


t+At 
=  min  /      J(x,0_,A)dA   +   1° 

6  t 


dx 

x   +  —  At,     t    +    At) 
—        at 


As   At   approaches   zero 


SI 


o      dx 


31 


I°(x,t)    =  min     J(x,6_,t)    At  +  J°(x,t)   +    (|j-)    dt   At  +  ft 


At 


1T°       '        dX  5T° 

i.e.       A*,e°.t>  +  (§-)    £  +  {§r-o 


"7 


which  may  be  written  as 


M  +  |f  =0 


(23) 


The  partial  differentiation  of  Equation  23  w.r.t.  x  yields 


Sx    3x-3t   — 


or 


3x     3t   — 


(24) 


However,  the  total  time  derivative  of  z  is 


dz 
dt 


oZ 


,  ,  dx 


.  j=  +  gq    - 

3t    ^x  ;   dt 


3J°    3(^'f)    ,3z«-  '  d^ 


dX 


)x      <-3x  j   dt 


(25) 


14 


ice 

If  +  3£  <£'D  -  0  (26) 

due  to  the  optimality  condition.   Expanding  Equation  25  gives  Equation 
19 ,  namely 

To    3f ' 


dz      3J 


so  that  partially  differentiating  w.r.t.  x  gives 


»2,    .2,      Ql   n     32f.    n     32f. 


(27) 


dt      3jc    82c  — 
Similarly , 

2 

3x      1  =  1       3x      *■     —  —      ; 


where 

se_* 

£  ■  -  te-  (29) 

and 

a2  n     32f. 

—   39_-  dx        39_  —    .£.   1  3_9_  3x_ 


To  evaluate  K  it  may  be  noted  from  Equation  17  that 


3J    3f* 

—  +  — —  z  =  0 

39    39   -   - 


^  ~2  +  ^TTsT  +  feH  A  zi~r+    I    zi  3iT^  +  3x~  tsT-J     =  ^ 


'-        3e-         —    -  "-      i=1      -    3Q-  i=i 


15 


i.e. 


30'   ,2T    n      ri 

^  3e2   i=i  1392 


and  K  = 


a2_    n     32f.  _1 

39    i=l     39 


(32) 


32J 
Therefore  it  is  possible  to  solve  for  — j  ,  namely 

3x 


„2T     dP    n     3*f.    r    3f  '     3f 
3x  i=l    3x     v    —       — 


(33) 


Now,  substituting  Equations  27  and  33  into  Equation  18  yields 


2T    n     3  f 


0    l       39    i=l    39 


69 


r  si  •     n     3f-,-  'i 

k   _     !=1    _    J 


+  69_'  R  <Sx  +  2"  6-'  -  -  &-\ 


(34) 


In  order  that  the  performance  index  will  converge  to  a  minimum,  the 
integral  in  Equation  34  must  be  less  than  zero,  i.e. 


16 


t  2 

few1  (H+  i  *i-T-)*±+  m  *  I  ^    «£ 

0         ^  86  i=l  80  <•      -  1=1      X     -       J 


+   60.'    •    R    •    6x  +  y  62L'    '    K    •    R    •    6x^   dt   <    0 


(35) 


In  addition,  the  convergence  idealy  should  be  as  fast  as  possible 
so  the  minimization  of  the  integral  is  considered: 


f  h  32T      n        9  f< 

(fix,  t)  -  /    I  66/  (li  +  J  «  _t)  59. 

t   ^    36   l-i   ae 


r  si  •       n  sf-  fi 


66 


+  69 '  R  6x  +  y  5x'  •  K  •  R  6x 


dA 


(36) 


Through  the  proper  choice  of  66_  and  denoting  the  minimum  by 

V  (6x,  t) ,  since  V(6x,  t)  as  given  by  Equation  36  is  quadratic  in  6x, 
the  minimum  of  V(6x,  t)  may  be  written  as  a  quadratic  expression,  as 


V  (6x,  t)  =  q(t)  +  (£(t))'  6x  +  6x'  Q(t)  6x 


(37) 


where    q(t)  =  scalar  function  of  t 

c^(t)  =  (nxl)  vector  function  of  t 

Q(t)  =  (nxm)  matrix  function  of  t  (symmetric) 

and   q(tf)  =  0    £(tf)  =  0     £(tf)  -  0  . 


From  Equation  37 


17 


3  V  (6x,  t)     d£(t)     dq_(t)  ' 


3t 


"  jl +  h )   <Sx  +  6xf  - 

dt      Mt   '        —         —  d 


dQ(t) 


dt 


6x 


(38) 


and 


3Vu(6x,  t) 
3x 


-  £(t)  +  2£(t)   6x 


(39) 


Minimization  of  V(6x,  t)  as  given  by  Equation  36  gives 


^•69°'  [3^J+  J 


z . 


ie2   i=i  x  3e2 


1 


56°  + 


—      1=1       —    J 


+  66°   R  •  6x  +  -|  6x »  K  R  6x 


dx 


+  £a1(t)  +  24x'a<t))«(53  +^iO 


dt'    dt 


d£(t)   '  d£(t) 

+  (- }   6x'  +  5x*  -: 6x  =  0 

vdx    '        —  —  dt     — 


(40) 


where 


2T    n     32f.  -1  t   ._ 
J  v  i\         \  c3J 


-(^\L^  ItM*^ 


3f' 


+  R  62L  +  far)  k(t)  +  ^(t)6*) 


■36 


(41) 


and 
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6dx 


1 


6f' 


f-(~)    «*+(^)  66.  . 


(42) 


When  the  optimal  control  as  given  by  Equation  41  is  substituted  into 
Equation  40  and  the  coefficients  of  6x   and  6x*    •    6x_   along  with  the  terms 
not  containing  6x   are  all  put  equal  to  zero  (to  satisfy  the  identify  for 
any  6x)  the  following  results: 


St^fs.^s  +  s.T-^lfwCtf)!-1* 


,  r3f  \  '   -1  r3f  S 


(43) 


d£       _  3f_'      3f_'  ' 

-R'  T^S  +  R'  T"1  (—h-  (— )   a 


dt 


3f '  '  3f '        3f ' 

!  fei-)    I"1  s  +  k=-)  t"1  (^)a 


I  81 


(44) 


and 


dQ 
dt 


3f 


3f '  »     3f_'       3f'  ' 


(45) 


where  S^  and  T  are  introduced  to  condense  the  notation  and  are  given  by 
n     3f 

I 

1-1 


(46) 


a2j       3  £ 


i-1     30 


(47) 


36    i-1 
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At  this  point  it  is  noted  that  the  matrix  (£  contributes  only  in- 
significantly to  the  control.   Furthermore  ^  appears  as  a  second  order 
term  itself.   Therefore,  to  facilitate  programming  on  the  digital 
computer  these  — = — -  equations  shall  be  discarded  and  ^  shall  be  put 

equal  to  zero.   Equation  43  is  not  required  for  the  evaluation  of  control; 
therefore  Equation  44  will  take  the  form 


|  -  R-  I"1  S  +  r   i1  (|f)  i-f)a  (48) 

and  the  change  in  the  control  simplifies  to 

-l  8i' 

6^=  -1       (l  +   R  62i  +  -39-  Sl)  •  <49) 


To  prevent  overstepping  in  control  adjustment,  Memiam  [23]  has 
suggested  the  introduction  of  a  constant  e  where  0  <  e  <_  1  in  Equation 
49  to  give 

66_  =  -eT_1[s  +  ||-  oj  -  T_1  R  6x  .  (50) 

Thus  it  is  now  worthwhile  to  detail  the  application  of  the  second  variation 
method  equations  developed  above. 

(1)  Assume  a  set  of  initial  value  for  9. 

(2)  Equations  1  and  2A  are  integrated  forward  from  t  =  0  to  t  =  tf; 
i.e.  (n+1)  equations  are  integrated  forward  in  time,  namely 

dx 


j£  =   J(x,  6,  t) 
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(3)  While  the  integration  is  carried  out,  the  values  of  x  are  retained 
in  the  computer  memory  at  small  time  intervals  to  approximate  the 
continuous  system. 

(4)  The  adjoint  Equation  19  plus  the  additional  Equations  28  and  48 
are  integrated  backwards,  i.e.,  2n+  ^^ — -  equations  are  integrated 
backwards  in  time  from  tf  to  0,  namely 


jz      jy_    (ill) 

dt     3x  "  ^3x  >    - 


dP     a2.    n     92f.    r  ,     3f»    ^ 

3x    i=l    3x     ^   —        —    ' 


^--T-i-'T-ffk-dfja. 


(5)  During  the  backward  integration,  the  values  of  T_,  S^,  £  and  R  are 
stored  in  computer  memory. 

(6)  The  new  value  of  control  is  calculated  from  Equation  50,  i.e. 

,a+«  .,«>  -  leT-^  +  g^^-Vv^x^  -x<J>) 

and  steps  2-6  are  carried  out  again. 

(7)  This  iteration  is  continued  until  no  further  change  in  6_  is  noticed 
or  until  the  performance  index  does  not  change.   The  former  is  more 
sensitive  [12]. 

If  the  performance  index  increases  during  some  iteration,  the  parameter 
e  is  halved  and  the  iteration  is  continued. 
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For  maximization  problems,  the  derivation  can  be  followed  on  the 
same  lines  and  it  will  be  seen  that  the  resulting  equations  are  the 


same. 
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2.4  Advantages  and  Disadvantages  of  the  Second  Variation  Method. 

The  foremost  advantage  of  the  second  variation  method  lies  in  its 
rapid  convergence.   Also,  unlike  the  first  variation,  the  optimum  can 
be  reached  with  reasonably  high  accuracy. 

The  theoretical  attractiveness  of  this  method,  however,  is  more  than 
offset  by  its  disadvantages.   First,  and  most  important,  the  initially 
assumed  trajectory  of  the  control  variable  must  be  sufficiently  close  to^ 
the  optimal  trajectory  for  convergence  to  be  obtained.   Second,  the 
number  of  equations  to  be  integrated  is  considerably  greater  than  required 
for  the  first  variation  method.   In  the  second  variation  method,  (n+1) 
equations  are  integrated  in  the  forward  direction,  i.e.  from  t=0  to  t=tf, 
and  (2n+n(n+l)/2)  equations  are  integrated  backwards  where  n  is  the  number 
of  state  variables  in  the  problem  under  consideration.   The  first  vari- 
ation method  requires  only  (n+1)  equations  to  be  integrated  in  the  forward 
direction  and  there  are  (n+1)  recursive  equations  in  the  backward  direction, 
Not  only  is  the  number  of  equations  involved  in  the  second  variation  method 
high  but  the  equations  themselves  are  more  complicated.   The  main  reason 
for  this  is  that  the  calculations  of  all  derivatives,  both  first  and  second 
order,  becomes  more  and  more  tedious  with  the  increasing  number  of  state 
and  control  variables.   All  the  multiplications  are  in  terms  of  matrices. 
Again,  the  inverse  of  T  has  to  be  computed  at  each  integration  step  in 
the  backward  integration.   Hence  the  programming  of  the  iteration  scheme 
with  the  required  equations  can  be  quite  complicated.   Instability  can 
arise  from  bad  starting  values,  i.e.  from  an  insufficiently  good  guess 
for  the  starting  trajectory  of  the  control  variable.   The  values  for  the 
parameter  e  have  to  be  established  by  trial  and  error  for  the  particular 
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problem.   The  higher  the  value  the  faster  the  convergence.   Finally,  this 
technique  cannot  handle  problems  involving  inequality  constraints. 
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3.   APPLICATIONS 

To  illustrate  the  use  of  the  second  variation  method,  three  numerical 
problems  in  the  field  of  production  planning  and  control  are  solved  in 
the  following  sections. 
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3.1  An  Inventory  Model 

The  Model 

The  following  is  a  simple  problem  in  the  field  of  production 
scheduling  and  inventory  control.   Assume  that  the  rate  of  sales  Q(t)  is 
known  with  certainty  and  that  the  rate  of  change  of  the  inventory  level 
I(t)  is  given  by 

« 

jfi&-   =  P(t)  -  Q(t)  (51) 

where  P(t)  is  the  production  rate  at  time  t.   The  problem  is  to  minimize 
the  cost  function. 

T 
CT  =  /  {CjC^  -  Kt))2  +  Cp  exp  (PM  -  P(t))2}  dt  (52) 

where  C   is  the  total  cost  of  inventory  and  production  and  C  is  the  min- 
imum production  cost  which  occurs  when  the  production  rate  equals  P  .   The 
quantity  P   can  be  considered  as  the  capacity  of  the  manufacturing  plant. 
Since  the  plant  is  designed  for  a  capacity  P  ,  an  increase  in  capacity 
may  require  additional  equipment  and  manpower  which,  due  to  contract 
agreements  cannot  be  reduced  easily.   I   can  be  considered  as  the  capacity 
for  the  storage  of  inventory  and  C   is  the  inventory  carrying  cost.   In 
many  practical  situations,  the  minimum  storage  cost  is  obtained  when  the 
storage  capacity  is  completely  filled.   Furthermore,  the  cost  function, 
Equation  (2) ,  has  the  smoothing  capability  which  is  frequently  desirable 
for  many  manufacturing  processes.   In  this  case,  I   and  P   can  be  con- 
sidered as  the  desirable  inventory  and  production  levels.   It  is  further 
assumed  that  the  sales  forecast  is  known  and  is  given  by  the  linear  re- 
lation 
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Q(t)  =  a+bt  (53) 

and  the  initial  inventory  is 

KO)  ■  c  (54) 

Recursive  Relations 

This  optimum  production  planning  problem  can  be  rewritten  into  the 
form  required  for  the  second  variation  method  as 
Let  xl(t)  =  I(t) 

6(t)  =  P(t)  . 

Equations  (51)  and  (54)  become 
dxl(t) 


dt 


=  6(t)  -  a  -  b(t)  (55) 


and  xl(0)  =  c  (56) 

Let 

x2(t)  =  /  CI(IM  -  HOT  +  Cp  exp  (PM  -  6(t))Z}dt  (57) 

Then 


x2(t)  =  CT  (58) 


^2U)  =  CjC^  -  xl(t))2  +  Cp  exp  (PM  -  G(t))2  (59) 

x2(0)  =  0  (60) 

Thus,  in  this  problem  there  is  one  state  variable,  namely  inventory  xl. 
The  control  variable  is  the  production  rate  9(t).   The  numerical  values 
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used  for   this  problem  are: 


a  =   2  b  =   l  c  =   5 


CT   -  0.1        I„  =   10  C„   =   0.001    . 

IMP 


PM  -   5  T  -   1 

M 


The  various  derivatives  required  for  obtaining  the  second  variational 
equations  are: 


£-lt=-2CI    ^M"Xl(t)) 


■&-  2  C 

ax2    * 


|f  =  -2  cp  exp  (pm-  e(t))2  •  (pm-  e(t)) 


^=2C  exp  (P  -  6(t))2   {1+2  (P  -  6(t))2} 
36 


^L-=  o 

36  3x 


3f!  32f- 

-  0 


3x  3x36 


32f 

71" ■  °  aT"  °     3^=  1 

3x  —  — 
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The  expression!*  for  the  terms  R,  s_,  T_  are: 

R  «=  P_  =  P         P_  being  1  dimensional 


8  =  -2  Cp  exp  (PM  -  0(t))2  (PM  -  G(t))  +  % 


T  =  2  Cp  exp  (pm  -  e(t))2  •  {1  +  2  (pm  -  e(t))2} 


The  second  variational  equations  (19,  28,  48)  become 


dF-f  =2CI  tln-x^t)]  (61) 


f£   -  £  =  -  2  Cz   +  P2  [2  Cp  exp  (PM  -  6(t))2  {1  +  2  (PM-6(t))2}] 


(62) 


d£F    dQF   P{-2C   exp(P   -  9(t))2(PM  -  6(t))  +  z  +  QF} 

JT  =  JT  = " f 5-    <63> 

2C   •  exp[P  -  e(t)T  {1  +  2[P  -6(t)]Z} 


and 


6<J+1>  .  e(J>  -  [E(.  +QF)](^  -  [P^Cx/^  -  ,  <*>»  . 


Thus  Equations  (61) ,  (62)  and  (63)  are  the  second  variational  equations 
and  Equation  (6A)  is  the  equation  for  finding  the  new  value  of  the  control. 


Table  1 


Effect  of  c  on  the  Rate  of  Convergence, 
of  Inventory,  0(t)  -  7,  x  (t)  -  5,  CKt<tf 
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Iteration 

E  =  0.1 

e  -  0.3 

e  -  0.5 

e  ■  0.7 

e  =  1.0 

1 

9.49995 

9.49995 

9.49995 

9.49995 

9.49995 

5 

9.38148 

9.39033 

9.37763 

9.36130 

9.34562 

10 

9.39130 

9.36100 

9.33515 

9.32687 

9.32586 

15 

9.38672 

9.33859 

9.32642 

9.32586 

9.32586 

20 

9.32649 

9.32585 

9.32588 

9.32586 

9.32586 

25 

9.32597 

9.32585 

9.32586 

9.32586 

30 

9.32587 

35 

9.32585 

40 

ii 

45 

ii 

50 

ti 

55 

it 

30 


10 
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Fig.  1   Convergence  >yate  of .Inventory. 
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Fig.  2.   Convergence  rate  of  Inventory. 
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Fig.  3  Convergence  rate  of  Inventory 


Table  1A 


Starting  Trajectories 
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0.00 

0.01 
0.02 
0.03 
0.04 
0.05 
0.06 
0.07 
0.08 
0.09 
0.10 
0.11 
0.12 
0.13 
0.14 
0.15 
0.16 
0.17 
0.18 
0.19 
0.20 
0.21 
0.22 
0.23 
0.24 
0.25 
0.26 
0.27 
0.28 
0.29 
0.30 
0.31 
0.32 
0.33 
0.34 
0.35 
0.36 
0.37 
0.38 
0.39 
0.40 
0.41 
0.42 
0.43 


0(t) 

7.189250000 
7.187477000 
7.079871000 
7.077939000 
7.104925000 
7.132057000 
7.149645000 
7.157608000 
7.158973000 
7.156805000 
7.153052000 
7.148717000 
7.144282000 
7.140039000 
7.135867000 
7.131478000 
7.127013000 
7.122532000 
7.118021000 
7.113436000 
7.108769000 
7.104050000 
7.099350000 
7.094633000 
7.089814000 
7.084888000 
7.079880000 
7.074810000 
7.069685000 
7.064513000 
7.059288000 
7.054013000 
7.048672000 
7.043266000 
7.037783000 
7.032225000 
7.026593000 
7.020878000 
7.015089000 
7.009217000 
7.003265000 
6.997230000 
6.991116000 
6.984918000 


x^t) 


5.000000000 
5.051840000 
5.103566000 
5.154113000 
5.204544000 
5.255143000 
5.305912000 
5.356759000 
5.407584000 
5.458322000 
5.508939000 
5.559419000 
5.609757000 
5.659949000 
5.709999000 
5.759905000 
5.809671000 
5.859291000 
5.908764000 
5.958094000 
6.007279000 
6.056316000 
6.105205000 
6.153948000 
6.202545000 
6.250991000 
6.299290000 
6.347438000 
6.395437000 
6.443284000 
6.490978000 
6.538519000 
6.585909000 
6.633147000 
6.680229000 
6.727155000 
6.773927000 
6.820544000 
6.867001000 
6.913302000 
6.959444000 
7.005426000 
7.051247000 
7.096908000 


Table  1A  (continued) 
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0.44 
0.45 
0.46 
0.47 
0.48 
0.49 
0.50 
0.51 
0.52 
0.53 
0.54 
0.55 
0.56 
0.57 
0.58 
0.59 
0.60 
0.61 
0.62 
0.63 
0.64 
0.65 
0.66 
0.67 
0.68 
0.69 
0.70 
0.71 
0.72 
0.73 
0.74 
0.75 
0.76 
0.77 
0.78 
0.79 
0.80 
0.81 
0.82 
0.83 
0.84 
0.85 
0.86 
0.87 
0.88 
0.89 
0.90 
0.91 
0.92 


6.9  78632000 
6.972250000 
6.965771000 
6.959211000 
6.952615000 
6.946091000 
6.939800000 
6.933926000 
6.928563000 
6.923484000 
6.917935000 
6.910558000 
6.899575000 
6.883111000 
6.859521000 
6.827771000 
6.787430000 
6.738802000 
6.682859000 
6.620663000 
6.553533000 
6.482723000 
6.409090000 
6.333541000 
6.256731000 
6.179157000 
6.101206000 
6.023172000 
5.945170000 
5.867554000 
5.790261000 
5.713470000 
5.637252000 
5.561565000 
5.486539000 
5.412208000 
5.338596000 
5.265617000 
5.193286000 
5.121719000 
5.050916000 
4.980895000 
4.911464000 
4.842929000 
4.775110000 
4.708012000 
4.641857000 
4.576363000 
4.511854000 


142408000 
187744000 
2329 J  4000 
27  7921000 
322764000 
367  440000 
411951000 
456297000 
500487000 
544522000 
7.588405000 
7.632135000 
7.675689000 
7.719036000 
7.762116000 
7.804861000 
7.847187000 
7.889011000 
7.930249000 
7.970827000 
8.010683000 
8.049766000 
8.088044000 
8.125484000 
8.162069000 
8.197786000 
8.232625000 
8.266588000 
8.299669000 
8.331870000 
8.363195000 
8.393647000 
8.423231000 
8.451951000 
8.479818000 
8.506833000 
8.533003000 
8.558340000 
8.582844000 
8.606528000 
8.629393000 
8.651453000 
8.672711000 
8.693175000 
8.712854000 
8.731755000 
8.749883000 
8.767251000 
8.783864000 
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Table  1A  (continued) 


0.93  4.448162000  8.799734000 

0.94  4.385418000  8.814865000 

0.95  4.323664000  8.829268000 

0.96  4.262953000  8.842953000 

0.97  4.203338000  8.855934000 

0.98  4.144891000  8.868217000 

0.99  4.087779000  8.879816000 

1.00  4.087779000  8.890743000 


Table  2 


Effect  of  e  on  the  Rate  of  Convergence 
of  Cost  Function  x  ,  0(t)  ■  7,  x  (t)  =-  5,  0<t<t- 
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Iteration 

c  =  0.1 

e  =  0.3 

e  =  0.5 

E  =  0.7 

e  =  1.0 

1 

0.95957 

0.95957 

0.95956 

0.95956 

0.95957 

5 

0.95232 

0.94536 

0.94392 

0.94356 

0.94342 

10 

0.94694 

0.94360 

0.94337 

0.94335 

0.94335 

15 

0.94498 

0.94339 

0.94335 

20 

0.94415 

0.94336 

M 

25 

0.94376 

0.94335 

it 

30 

0.94356 

ii 

n 

35 

0.94347 

it 

ii 

40 

0.94342 

ii 

M 

45 

0.94339 

ii 

n 

50 

0.94339 

ii 

M 

n 

55 

0.94336 

ii 

n 

ii 
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Fig.  4  Convergence  rate  of  Cost  Function. 
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Fig.  6   Convergence  rate  of  Cost  Function, 


Table  3 


Effect  of  e  on  Rate  of  Convergence  of  the 
Production  Rate  8,  8(t)  -  7,  x  (t)  ■  5,  CKt<_t 


40 


Iteration     c  =  0.1 

e  =  0.3 

e  =  0.5 

c  =  0.7 

E  =  1.0 

1 

B(t03 

7.03101 

7.09304 

7.15507 

7.21709 

7.31013 

e(tf) 

6.97778 

6.93333 

6.88889 

6.84444 

6.77778 

5 

e(t03 

7.10105 

7.17094 

7.18717 

7.18952 

7.18934  • 

e(tf: 

6.88690 

6.64770 

6.38816 

6.10440 

5.62812 

10 

e(tQ: 

7.14225 

7.18636 

7.18927 

7.18933 

7.18933 

e(tf: 

6.76850 

6.23612 

5.58561 

5.06248 

5.0000 

15 

8(t0] 

7.16293 

7.18884 

7.18933 

7.18933 

7.18933 

e(tf) 

6.64410 

5.74570 

5.03792 

5.00015 

5.00000 

20 

e(tQ) 

7.17416 

7.18925 

7.18933 

7.18933 

7.18933 

6(tf) 

6.51294 

5.25530 

5.00119 

5.00000 

5.00000 

25 

e(t0: 

7.18051 

7.18932 

7.18933 

7.18933 

7.18933 

e(tf) 

6.37410 

5.04754 

5.00004 

5.00000 

5.00000 

30 

8(t0] 

7.18416 

7.18933 

7.18933 

7.18933 

7.18933 

e(tf: 

6.22667 

5.00802 

5.00000 

5.00000 

5.00000 

35 

e(t0: 

7.18629 

7.18933 

7.18932 

7.18933 

7.18933 

e(tf: 

6.06984 

5.00135 

5.00000 

5.00000 

5.00000 

40 

e(t0: 

7.18754 

7.18933 

7.18933 

7.18933 

7.18933 

e(tf; 

I     5.90348 

5.00023 

5.00000 

5.00000 

5.00000 

45 

e(t0: 

7.18828 

7.18932 

7.18933 

7.18933 

7.18933 

e(tf; 

1     5.72933 

5.00003 

5.00000 

5.00000 

5.00000 

50 

8(t0; 

)      7.18871 

7.18932 

7.18933 

7.18933 

7.18933 

e(tf" 

)      5.55353 

5.00000 

5.00000 

5.00000 

5.00000 

55 

6(to 

7.18896 

7.18932 

7.18933 

7.18933 

7.18933 

6(tf 

>     5.38931 

5.00000 

5.00000 

5.00000 

5.00000 
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Fig,  7  Convergence  rate  of  Production  Rate. 
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Fig.  8   Convergence  rate  of  Production  Rate. 
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Fig.  9   Convergence  rate  of  Production  rate. 
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NufflPI   I  Cl\\  lllH 

This  problem  was  solved  by  two  approaches.   In  the  first,  the 
second  variation  method  was  used  in  combination  with  the  first  variation 
method.   The  same  problem  is  solved  by  Lee  and  Shaikh  (20).   The  values 
of  the  state  and  control  variables  (at  all  grid  points)  were  taken  from 
the  results  of  the  first  variation.   In  particular,  the  values  of  Xl(t) 
and  6(t)  were  taken  from  the  21st  iteration  of  the  first  variation  and 
fed  as  good  starting  values  for  the  second  variation.   These  values  are 
listed  in  Table  1A.   In  the  second  approach,  the  second  variation  was 
tried  directly  by  itself.   For  this,  a  guess  was  made  for  the  starting 
values  of  the  state  variable  and  the  control  variables. 

An  interesting  parameter  in  the  computation  is  the  step  size  e  which 
determines  the  magnitude  of  the  step  taken  in  each  iteration.   In  the 
solution  of  this  problem,  a  series  of  values  of  e  were  selected  and  the 
computation  was  carried  out  for  each.   Tables  1,  2  and  3  show  the  con- 
vergence rate  of  inventory  xl,  cost  function  x2  and  the  production  rate 
6(t),  respectively.   These  tables  are  for  a  constant  starting  value  of 
the  control  variable,  namely  6(t)  =  7,  0  <_  t  <  t.  and  a  constant  starting 
value  of  the  state  variable,  namely  xl(t)  =  5,  0  <_  t  <_   tf .   It  is  seen 
that  for  e  =  1,  the  fastest  convergence  rate  is  obtained  while  the  con- 
vergence rate  slows  down  when  its  value  is  decreased.   Figures  1  through 
9  show  the  rate  of  convergence  of  the  inventory,  production  rate  and  the 
cost  function  for  different  values  of  e  . 

Regarding  the  starting  trajectory  of  the  control  variable,  it  was 
found  that  only  the  constant  trajectories  between  6(t)  ■  7  and  6(t)  =  8 
would  lead  to  convergence.   For  all  other  control  variable  values,  the 
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the  problem  would  not  converge.   These  other  values  were: 

0(t)  ■  1,2,3,4,5,6    and  0(t)  =  9,10,11,    0_Stltf 

Also,  the  combination  of  the  first  and  the  second  variation  required  about 
50  iterations  to  reach  the  optimal  with  e  ■  0.3.   A  higher  value  of  e  could 
not  be  used  as  it  led  to  overstepping  in  this  situation. 
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3.2   An  Inventory  and  Advertising  Model 
The  Model 

This  model  is  an  extension  of  the  one  formulated  by  Teichroew  (27). 
Consider  a  marketing  situation  where  only  a  certain  number  of  possible 
customers  possess  certain  information  about  a  firm's  product.   Suppose 
that  the  total  number  of  such  possible  customers  remains  constant  and 
that  the  diffusion  of  information  occurs  only  through  personal  contact. 
The  number  of  contacts  made  by  an  informed  person  in  a  unit  time  is  known 
as  contact  coefficient.   In  a  contact,  the  contactee  receives  information 
if  he  does  not  already  have  it;  if  he  already  has  it,  the  contact  is 
wasted  insofar  as  increasing  the  number  of  informed  persons  is  concerned. 

Let  K(0)  =  K„  =  number  of  informed  persons  at  time  tn 

N  =  total  number  of  persons 

c  =  contact  coefficient,  the  number  of  contacts  made  by 
one  informed  person  per  unit  time 
K(t)=  number  of  informed  persons  at  time  t. 
Then  K(t)/N  =  proportion  of  informed  persons  at  time  t 

1  -  K(t)/N  =  proportion  of  uninformed  persons  at  time  t 
c.K(t).dt  =  contacts  made  during  a  time  interval  dt. 

Clearly  dK(t)  =  c.K(t) .dt . (l-K(t) /N) 

Thus  the  equation  governing  the  process  is 

^^-     =  c.K(t).(l-K(t)/N)  (65) 

at 

Suppose  next  that  the  firm  can  influence  the  number  of  contacts  by 

spending  money  on  advertising.   In  particular  it  can  increase  the  number 
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of  contacts  made  by  the  informed  persons  (above  the  ones  included  in  c) 
by  an  additional  number  A  per  unit  time. 
Equation  (65)  now  becomes 

|p-}-  =  K(t).(c+A(t)).(l-K(t)/N)  (66) 

If  each  successful  contact  results  in  the  sale  of  n  units  of  the 
firm's  product  and  if  Q(t)  represents  the  sale  at  time  t,  then 

Q(t)  =  n  K(t) 
Letting  n=l  and  substituting  Q(t)  for  K(t)  in  Equation  (66),  then 

f^  =  Q(t).(c+A(t)).(l-Q(t)/N)  (67) 

The  rate  of  change  of  the  firm's  inventory  is  given  by 
dX(t) 


dt 


=  P(t)  -  Q(t)  (68) 


where  P(t)  =  production  rate  at  time  t. 

The  production  rate  is  assumed  to  be  a  linear  function  of  time 

P(t)  =  a+bt  (69) 

where  a  and  b  are  constants. 

This  assumption  is  made  to  simplify  the  model  by  avoiding  a  second 
control  variable. 

The  firm's  management  wishes  to  maximize  the  profit 

T 
ST  =  /  [F.Q(t)  -  Cx  (Px  -  x(t))2  -  CA  A2(t)  Q(t)]dt  (70) 

where  S   is  the  total  net  profit. 
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F  is  the  revenue  from  Lhe  .sale  of  one  unit  of  the  product.   C   is 
the  inventory  carrying  cost  and  has  the  same  significance  as  in  the 
model  descrihed  in  Section  3.1.   P_  can  be  considered  as  the  capacity  for 
the  storage  of  inventory.   C   is  the  cost  of  advertising. 

Equations  (67)  through  (70)  represent  the  system  under  consideration. 
The  system  has  two  state  variables,  inventory  X(t)  and  sales  Q(t),  and 
there  is  one  control  variable,  advertising  A(t). 

* 

The  initial  conditions  and  the  numerical  values  used  are: 

a  =  0.7   b  =  1.0   c  =  2.0   N  =  1.5   F  =  10.0 

C  =  0.15   P  -  1.0   CA  =  1.0   X(0)  =  0.2   Q(0)  =  0.2 

Recursive  Relation 

The  necessary  relations  for  the  second  variation  can  be  obtained  in 
the  following  manner.   Note  that  in  these  derivations  x(t)  denotes  the 
state  variable  vector  while  x(t)  denotes  the  inventory.   From  Equation  70, 
then, 


J  =  Q.F  -  CI(P].  -  x(t))2  -  CAQA2(t). 


The  various  derivatives  required  for  obtaining  the  second  variation  equations 
are: 


r 


3J 

ox 


*1 

3x 


3Q 


2CI(PI  -  x(t)) 


F  -  CAA*(t) 
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3*2 


_3_ 
3x 


3  J 


3x 


3_ 
3Q 


fajl 


.  3x, 

V.  J 


3 
3x 

fajl 

3_ 
3Q 


3J 


[3Q 


-2C  0 


If  -  set  ■  "2ca  «(c)  A(t> 


^f  ■=  -^-  -  -2  c    q(t) 

38"'        3A   (t)  A 


32J 


36    3x        3A(t) 


3J 

3x 


M 
3Q 


-    [0 


-2CA  A(t)] 


3x 


0 


ch 


3x    36 


=      [0      ,        0] 


3f. 
3x~ 


C+A(t) 


1  - 


2Q(t) 

N 


32£. 


3x    36_ 


o,    (i-^I) 
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30~ 


-   0 


DO2 


3f. 

Tq 


Q(t) 


i  _aui 

N 


S2£. 


39' 


Ml 
3x 


3x    L    1J 


—    [f    1 
3Q    L    1J 


—  [f    ] 
3x    l    2J 


-2-   ff    1 
3Q    L    2J 


-1 


[C+A(t)]jl-^ 


32f. 


3x 


-2f 


3x 


^  (C  +  A(t)) 


3f  ' 
36 


3A(t) 


[f. 


f2] 


Q(0 


1  - 


Q(t) 

N 


Expressions  for  the  terms  R,  s_,  T_   from  Equation  30  result  in 
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R  =   [0     -2CA  A(t)  ]  + 


0      Q(t) 


N 


P        P 
rll       12 


P        P 
12       22 


z2(l 


2Q(t) 


[0 


■2CA  A(t)  ] 


P12  Q(t)(l 


Q(t) 

N 


)  , 


P22  Q(t)(l 


Q(t) 

N 


z2  (1 


2Q(t) 

N 


Let    R  =  [R1,    R2] 


where 


\   =  P12   Q(t) 


1  - 


Q(t) 

N 


R2  =  "2  CA  A(t)  +  P22  Q(t) 


i  .asa. 

N 


+  z. 


1  - 


2Q(t) 


Equation  46  gives 


S  =  -2CA  Q(t)  A(t)  +  z2  Q(t) 


N 


and  Equation  47  gives, 


T  -  -2CA  Q(t). 

It  is  now  possible  to  determine  the  2n  +  — i.e.  (2+2+3)  or 

seven  equations  to  be  integrated  backwards.   Equation  (19)  becomes 
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dt 


-2CX  (Pj  -  x(t)) 


-F  +  CA  A  (t) 


-1 


[C+A(t)] 


\        2Q(t)^ 

N 


1 

rv 

i 

Z2 

-2^  (Pj  -  x(t) 


-F  +  C  A  (t) 


-Zx  +  z2  [C  +  A(t)] 


1  - 


2Q(t) 

N 


dz 

—  =-2Cl  [Px  -x(t)] 


(71) 


dz 
d 


^  -  -F  +  CA  A2(t)  +  z±   -   z2  [C  +  A(t)]|l  -  ^9ili 


(72) 


Thus  Equations  (71)  and  (72)  correspond  to  Equation  (19).   Equation  (28) 
in  this  case  becomes 
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2C, 


(    0 


dP 
dt 


2z, 

~n" 


[C  +  A(t)] 


Pll  +  P12  [C  +  A(t)] 


1  -Mil 

N 


Pll  +  P12  [C  +  A(t)] 


1  - 


2Q(t) 

N 


-  2P12  +  2P22  [C  +  A(t)] 


1  .  Mil 

N 


Hence  Equation  (28)  is  represented  by  the  following  three  equations: 


dPll  2 

—  "  2CI  +  R1T 


(73) 


dP 


12 


-P.-  -  P10  [C  +  A(t)] 


dt      11    12 


1-^1  +  RlR2T 

I       N   i     1   2 


(7A) 


dP22    2. 


dt 


■jj—  [C  +  A(t)]  +  2P12 


2P22  [C  +  A(t) 


2Q(t)) 

N 


+  R2T 


(75) 


To  avoid  confusion,  the  £  in  the  derivation  of  the  method  given  in  Equation 
(48)  is  denoted  by  QF  here.   Thus  Q  still  represents  the  sales  for  this 


problem. 

Equation  (48)  is  given  by 
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42, 
dt 


I 
Rl 


V      J 


(f)  s  + 


f        \ 

/ 

h 

(|) 

lR2, 

V 

Q(t) 


.  2iO 


' 

QF, 

. 

,QF2. 

-1 


[C  +  A(t)] 


1  - 


2Q(t) 


QF 


QF, 


f         1 

f 

R  s 

H 

- 

R2s 

T 
*> 

. 

Rl 

■~   •  QF2  •  Q(t) 


N 


QF, 


•  QCt)  -  (l  -  *tf 


I 


-  QF, 


QF2  [C  +  A(t)] 


1  .  MEi 

N 


dQF.    R,s    R 

-  -=-  +  —  •  QF„  •  Q(t) 

dt      T     T   w  2   wv  ' 


1  - 


Q(t) 


+  QF, 


(76) 
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Table  4 


Starting  Trajectories 


C  Kt)  Q(t) 

0.00  0.200  0.200 

0.01  0.203  0.200 

0.02  0.205  0.201 

0.03  .  0.208  0.208 

0.04  0.210  0.212 

0.05  0.211  0.220 

0.06  0.212  0.228 

0.07  0.217  0.230 

0.08  0.220  0.238 

0.09  0.222  0.245 

0.10  0.226  0.250 

0.11  0.228  0.261 

0.12  0.230  0.272 

0.13  0.233  0.283 

0.14  0.238  0.300 

0.15  0.242  0.317 

0.16  0.250  0.339 

0.17  0.252  0.410 

0.18  0.260  0.430 

0.19  0.264  0.450 

0.20  0.270  0.460 

0.21  0.274  0.483 

0.22  0.280  0.503 

0.23  0.284  0.520 

0.24  0.290  0.540 

0.25  0.293  0.560 

0.26  0.300  0.580 

0.27  0.301  0.600 

0.28  0.306  0.620 

0.29  0.310  0.648 

0.30  0.318  0.665 

0.31  0.320  0.690 

0.32  0.324  0.702 

0.33  0.330  0.724 

0.34  0.336  0.745 

0.35  0.340  0.760 

0.36  0.342  0.785 

0.37  0.346  0.800 

0.38  0.350  0.814 

0.39  0.351  0.830 

0.40  0.355  0.842 

0.41  0.359  0.857 

0.42  0.360  0.870 

0.43  0.362  0.880 

0.44  0.368  0.890 
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0.45 


Table   4      (continued) 


0.46  0.371 

0.47  0.372 


0.370  0.900 

0.910 


°'54  0.389 

°'55  0.390 

0.56  0.391 

0.57  0.392 


0.72 

0.73  0.422 

°-74  0.426 

°-75  0.430 

°'76  0.435 

°-77  0.438 

°-78  0.440 

°-79  0.441 

0.80  0.447 

°-81  0.450 

°-82  0.451 

°-83  0.453 

°-84  0.460 

°-85  0.461 

°-86  0.462 

°-87  0.468 

°-88  0.470 

°-89  0.473 

°-90  0.475 

°-91  0.480 


0.92 

°-93  0.488 


0.920 
0.923 


o'f  °-373 

o'ln  °'378  0.930 

0*5?  °'380  0.936 

Q\?  °-381  0.940 

0*53  n'382  0-945 

°-3  0.950 

0.951 
0.958 
0.962 
0.970 


o'll  °*393                                                            0.971 

0.398  n  q7o 

°-60  o  Ann 

n  A1  °-400                                                    0.982 

olo  °-40°                                                              0-990 

?  O'401                                                                        0    99S 

0-63  n  Am 

n  (,l  3                                                    i-OOO 

0*65  O'408                                                               1.002 

0-65  0.410 

0.66  0.411 


1.010 
1.015 


0-412  1.020 


1.028 


0.67 

0-68  0.416 

0-69  0.417 

°-7°  0.419  i*03? 

0-421  '  1.045 

1.050 


1.051 

1.056 

1.060 

1.065 

1.070 

1.071 

1.078 

1.081 

1.090 

1.092 

1.100 

1.101 

1.110 

1.112 

1.120 

1.130 

1.135 

1.140 


0.482  lel43 

1.150 
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Table  4   (continued) 


0.94  0.490  1.160 

0.95  0.496  1.165 

0.96  0.500  1.170 

0.97  0.503  1.173 

0.98  0.510  1.180 

0.99  0.513  1.190 

1.00  0.519  1.200 
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dt      T    T    Q  2 


Q(t) 


1  - 


Q(0 


-  QF2  [C  4-  A(t)]  • 


1  _  Mil 


(77) 


Thus  Equations  (76)  and  (77),  represent  Equation  (48).   The  equation  for 
improving  the  control  variable  becomes 


A(t) 


(J+1)-A(t)^+f 


s  +  QF_  •  Q(t) 


1  .act* 


i{R1(x(J+1).x^)+R2(Q<J+1>.Q^>)} 


Equation  (78)  represents  Equation  (50)  (78) 

This  problem  illustrates  how  tedious  the  calculations  become  when  the 
number  of  variables  increases. 

Numerical  Results 

In  here,  the  starting  trajectories  of  the  two  state  variables, 
inventory  I(t)  and  the  sales  Q(t),  were  fed  from  the  results  of  the  solution 
of  the  same  problem  by  dynamic  programming.   These  values  are  listed  in 
Table  4.   Actually  these  values  are  obtained  after  dividing  the  original 
results  by  100.   This  was  required  to  prevent  the  exponential  overflow  of 
the  system  of  equations.   The  starting  trajectory  of  the  control  variable 
was  tried  in  the  range  of  0.001  to  6.0.   It  was  found  that  all  these  values 


would  work;  however,  the  best  value  was  found  to  be  0(t)  -  0.5,  0<t<t.a 
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Table  5 


Effect  of  e  on  the  Rate  of  Convergence 
of  I(tf)  with  AQ(t)  -  0.5. 


Iteration 

E  =  0.1 

E  =  0.3 

e  =  0.5 

E  =  0.7 

1 

0.8524 

0.8524 

0.8524 

0.8524  . 

5 

0.7137 

0.6274 

0.6076 

0.6277 

10 

0.6546 

0.5990 

0.5939 

0.5929 

14 

0.6307 

0.5948 

0.5935 

0.5934 

16 

0.6227 

0.5941 

M 

0.5935 

17 

0.6194 

0.5939 

H 

ii 

21 

0.6096 

0.5936 

ii 

ii 

The  Values  of  IQ(t)  &  QQ (t)  are  obtained  from  Table  4. 
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Fig.  10    Convergence  rate  of  Inventory. 
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Fig.  11   Convergence  rate  of  Inventory, 
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Fig.  12   Convergence  rate  of  Inventory. 
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Table  6 


Effect  of  g  on  the  Rate  of  Convergence 


of  Q(t  ) ,  with  A  (t)  =  0.5. 


Iteration 

c   =  0.1 

e   =  0.3 

e   =  0.5 

e   =  0.7 

1 

0.9781 

0.9781 

0.9781 

0.9781 

5 

1.135 

1.198 

1.206 

1.172 

10 

1.179 

1.218 

1.222 

1.223 

14 

1.197 

1.221 

ii 

ii 

16 

1.202 

1.222 

ii 

ii 

17 

1.205 

n 

ii 

ii 

21 

1.211 

ti 

ii 

ii 

The  Values  of  IQ(t)  &  Q0(t)  are  obtained  from  Table  4. 
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Fig.  13   Convergence  rate  of  Sales. 
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Table  7 


Effect  of  e  on  the  Rate  of  Convergence 
of  Total  Profit,  with  AQ(t)  -  0.5. 


1 teration 

E     =    0.1 

e   -   0.3 

e   =  0.5 

e    =   0.7 

1 

5.298 

5.298 

5.298 

5.298       . 

5 

6.260 

6.596 

6.621 

6.571 

10 

6.527 

6.626 

6.626 

6.626 

14 

6.589 

ti 

ii 

ii 

16 

6.604 

ii 

ii 

ii 

17 

6.609 

it 

ii 

M 

21 

6.620 

ii 

ii 

ii 

The  Values  of  IQ(t)  &  QQ(t)  are  obtained  from  Table  4. 
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Fig.  16   Convergence  rate  of  Profit 
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Fig.  17  Convergence  rate  of  Profit, 
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Fig.  18   Convergence  rate  of  Profit. 
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Table  8 


Effect  of  e  on  the  Rate  of  Convergence 
of  A(t  ) ,  with  A  (t)  =  0.5. 


Iteration 

e  =  .01 

e  =  0.3 

e   =  0.5 

e  -  0.7 

1 

1.320 

2.960 

4.599 

6.239 

5 

3.088 

4.841 

5.218 

5.269 

10 

4.091 

5.174 

5.222 

5.221 

14 

4.525 

5.213 

5.221 

ii 

16 

4.672 

5.217 

ii 

ii 

17 

4.733 

5.219 

M 

ii 

21 

4.917 

5.220 

ii 

ii 

The  Values  of  IQ(t)  &  QQ(t)  are  obtained  from  Table  4. 
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Fig.  19   Convergence  rate  of  Advertisement. 
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Fig.  20  Convergence  rate  of  Advertisement. 
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Fig,  21   Convergence  rate  of  Advertisement, 
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Table  9 


Starting  Trajectories  for  Inventory,  Sales 
and  Advertisement,  0^t<_tf  . 


I0(t)  -  0.2  2     IQ(t)  «=  0.5 

Q0(t)  -  0.2  QQ(t)  =  0.5 

AQ(t)  =  0.5  AQ(t)  =  2.0 


IQ(t)  =  0.5  4     IQ(t)  =  0.6 

Q0(t)  =  1.0  QQ(t)  =  1.3 

AQ(t)  =  2.0  AQ(t)  =  2.0 


IQ(t)  =  0.6 
Q0(t)  =  1.3 
AQ(t)  =  5.0 
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Table  10 


Effect  of  e  on  Rate  of  Convergence  of 
I(tf)  with  IQ(t)  =  QQ(t)  =  0.2,  AQ(t)  =  0.5,  0<t<tf< 


Iterat: 

.on 

c  =  0.1 

c  =  0.3 

e  =  0.5 

e  =  0.7 

1 

0.852U 

0.852U 

0.852U 

0.852U 

5 

0.726U 

0.63^3 

0.611U 

0.6322 

10 

0.662U 

0.5999 

0.59^0 

0.5928 

15         0.6309       0.59^5        0.5935        0.5935 
20         O.61U2       0.5936 
25         0.6051       0.5935 
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Fig.  22   Convergence  rate  of  Inventory. 
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Fig.  23  Convergence  rate  of  Inventory. 


Table  11 


Effect  of  e  on  Rate  of  Convergence  of 
Q(tf),  IQ(t)  =  QQ(t)  -  0.2,  AQ(t)  -  0.5,  0^t£tf. 
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Iteration 

e  =  0.1 

e  =  0.3 

£  =  0.5 

e  «  0.7 

1 

0.9781 

0.9781 

0.9781 

0.9781  . 

5 

1.083 

1.178 

1.198 

1.165 

10 

1.153 

1.215 

1.222 

1.223 

15 

1.185 

1.221 

ii 

1.222 

20 

1.202 

1.222 

ii 

ii 

25 

1.211 

ii 

ii 

ii 
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Fig.  24   Convergence  rate  of  Sales. 
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Fig.  25  Convergence  rate  of  Sales. 
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Table  12 


Effect  of  e   on  Rate  of  Convergence  of 
Total  Profit,  IQ(t)  -  Q  (t)  «=  0.2,  AQ(t)  -  0.5,  0<t<tf. 


Iteration 

e  =  o.i 

E  =  0.3 

e  =  0.5 

e  =  0.7 

1 

5.298 

5.298 

5.298 

5.298  . 

5 

6.246 

6.586 

6.614 

6.554 

10 

6.518 

6.626 

6.626 

6.626 

15 

6.595 

it 

ii 

ii 

20 

6.617 

M 

H 

ii 

25 

6.624 

ii 

ii 

n 
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Fig,  26   Convergence  rate  of  Profit 
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Fig.  27    Convergence  rate  of  Profit. 


Table  13 


Effect  of  e  on  Rate  of  Convergence  of 
A(tQ),  IQ(t)  ■  QQ(t)  =  0.2,  AQ(t)  =  0.5,  0<t<tf 
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Iteration 

E  =  0.1 

E  =  0.3 

e  =  0.5 

e  =  0.7 

1 

1.320 

2.960 

4.599 

6.239  . 

5 

2.973 

4.804 

5.223 

5.284 

10 

4.005 

5.169 

5.223 

5.220 

15 

4.549 

5.215 

5.221 

5.221 

20 

4.847 

5.220 

it 

it 

25 

5.012 

5.221 

it 

ti 
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Pig,  28   Convergence  rate  of  Advertisement. 
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Fig.  29   Convergence  rate  of  Advertisement. 
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Table  14 


Effect  of  Different  Starting  Trajectories 
on  I(tf)  with  E  =  0.4. 


IQ(t)  =  0.5    IQ(t)  =  0.5    IQ(t)  =  0.6    IQ(t)  =  0.6 

Iteration    QQ(t)  =  0.5    QQ(t)  =  1.0    QQ(t)  =  1.3    QQ(t)  =  1.3 

AQ(t)  =  2.0    AQ(t)  =  2.0    AQ(t)  =  2.0    AQ(t)  =  5.0 


1 

0.6134 

0.6134 

0.6134 

0.3305 

5 

0.6024 

0.5922 

0.5876 

0.5356 

10 

0.5945 

0.5939 

0.5936 

0.5887 

15 

0.5936 

0.5936 

0.5935 

0.5932 

20 

0.5935 

0.5935 

ii 

0.5934 

Table  15 

Effect  of  Different  Starting  Trajectories  on 
Rate  of  Convergence  of  Q(tf)  with  e  =  0.4. 
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Iteration 


IQ(t)  =  0.5 
Q0(t)  =  0.5 
AQ(t)  =  2.0 


I0(t)  =  0.5 
Q0(t)  =  1.0 
AQ(t)  =  2.0 


IQ(t)  =  0.6 
Q0(t)  =  1.3 
AQ(t)  =  2.0 


IQ(t)  =  0.6 
Q0(t)  =  1.3 
AQ(t)  =  5.0 


1 

1.340 

1.340 

1.340 

1.491 

5 

1.219 

1.243 

1.255 

1.316 

10 

1.222 

1.224 

1.225 

1.232 

15 

it 

1.222 

1.222 

1.223 

20 

it 

ii 

it 

1.222 
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Table  16 

Effect  of  Different  Starting  Trajectories  on 
Rate  of  Convergence  of  Total  Profit  with  e  =  0.4, 


IQ(t)  =  0.5    IQ(t)  =  0.5    IQ(t)  =  0.6    IQ(t)  =  0.6 

Iteration    QQ(t)  =  0.5    QQ(t)  =  1.0    QQ(t)  =  1.3    QQ(t)  =  1.3 

AQ(t)  =  2.0    AQ(t)  =  2.0    AQ(t)  =  2.0    AQ(t)  =  5.0 


1 

4.668 

4.668 

4.668 

-16.120 

5 

6.621 

6.588 

6.551 

6.249 

10 

6.626 

6.625 

6.625 

6.621 

15 

it 

6.626 

6.626 

6.626 

20 

ii 

M 

ii 

it 
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Table  17 


Effect  of  Different  Starting  Trajectories  on 
Rate  of  Convergence  of  A(t~)  with  e  =  0.4. 


IQ(t)  =  0.5    IQ(t)  =  0.5    IQ(t)  =  0.6    IQ(t)  =  0.6 

Iteration    .QQ(t)  =  0.5    QQ(t)  =  1.0    QQ(t)  =  1.3    QQ(t)  =  1.3 

AQ(t)  =  2.0    AQ(t)  =  2.0    AQ(t)  =  2.0    AQ(t)  =  5.0 


1 

2.330 

1.442 

-1.176 

1.823 

5 

4.910 

4.677 

4.556 

4.476 

10 

5.210 

5.182 

5.173 

5.151 

15 

5.220 

5.218 

5.217 

5.215 

20 

5.221 

5.221 

5.221 

5.220 
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Fig.  30   Convergence  rate  of  Advertisement, 
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Fig.  31   Convergence  rate  of  Advertisement. 
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Fig.  32   Convergence  rate  of  Advertisement. 
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The  results  Starting  with  this  trajectory  were  explored  in  detail  with 
different  value  of  c  .   Tables  5  through  8  show  the  convergence  rate  of 
inventory,  sales,  profit  function  and  advertisement,  respectively,  for  the 
different  values  of  c. 

Figures  10  through  12  show  the  convergence  rate  of  inventory  for  dif- 
ferent values  of e  .  Similarly,  Figs. 13  through  15  show  the  convergence 
rate  of  sales,  Figs.  16  through  18  of  profit  function  and  Figs.  19  through 
21  of  advertisement  for  different  values  of  e  .  The  maximum  value  of  e  that 
would  lead  to  convergence  in  this  case  was  found  to  be  0.7.  c  =1.0  would 
lead  to  exponential  overflow  in  this  situation.  Another  interesting  point 
noted  was  that  almost  the  same  convergence  rate  was  obtained  with  e  =  0.5 
and  with  e  =  0.7.   Thus  a  higher  £  did  not  increase  the  convergence  rate. 

In  an  another  approach  to  this  problem,  a  number  of  different  starting 
trajectories  for  inventory,  sales  and  advertisement  were  used.   These  are 
listed  in  Table  9.   Set  (1)  of  the  various  trajectories  listed  in  Table  9 
was  explored  in  detail  with  different  values  of  c  . 

Tables  10  through  13  show  the  convergence  rate  of  inventory,  sales, 
profit  function  and  advertisement  respectively,  for  different  values  of z  . 
Figures  22  and  23  show  the  convergence  rate  of  inventory  for  different 
values  of  c  .   Similarly  Figs.  24  and  25  show  the  convergence  rate  of 
sales,  Figs.  26  and  27  of  profit  function  and  Figs.  28  and  29  of  advertise- 
ment for  different  values  of  e  . 

The  remaining  starting  trajectories  from  Table  9,  namely  sets  (1) 
through  (5)  were  tried  with  e  =  0.4.   Tables  14  through  17  list  the  con- 
vergence rate  of  advertisement  for  these  trajectories.   Figs.  30  through 
33  show  the  convergence  rate  of  advertisement  for  these  different  tra- 
jectories. 
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The  starting  trajectories  (1)  through  (5)  from  Table  9  led  to  con- 
vergence almost  in  the  same  number  of  iterations.  The  maximum  value  of 
e  that  would  lead  to  convergence  was  found  to  be  0.7  in  this  case  also. 

Thus  it  is  seen  that  the  problem  is  very  stable  and  that  the  optimum 
can  be  reached  almost  with  any  reasonable  values  of  the  starting  tra- 
jectories. 

Another  computational  feature  that  was  encountered  in  the  solution 
of  this  problem  was  regarding  the  numerical  solution  of  the  differential' 
equations.   As  their  number  increased  it  was  found  advisible  to  use  the 
IBM  subroutine  "RKGS"  for  their  numerical  solution.   However,  this  sub- 
routine imposed  the  problem  of  accuracy  which  has  to  be  specified  by  the 
user.   This  is  the  accuracy  against  which  the  results  are  checked  after 
each  integration  step.   If  the  accuracy  is  too  low,  the  integration  step 
size  is  halved  and  this  continues  until  the  specified  accuracy  is  obtained. 
Thus,  if  the  accuracy  is  not  appropriate,  the  grid  points  may  not  be  the 
ones  desired  by  the  user.   The  calculations  of  R,S_,    and  T_  should  be  done 
both  in  subroutine  "FCT"  and  "OUTP".   (See  appendix  7.2)   Also,  to  test 
the  fact  that  this  method  would  lead  to  convergence  at  the  nearest  stationary 
point  regardless  of  whether  it  is  a  maximization  or  a  minimization  problem, 
the  objective  function  was  made  negative  and  the  same  problem  solved  again. 
The  results  agree  in  both  the  cases.   Thus  whether  a  maximum  or  minimum 
will  be  reached  all  depends  on  the  nature  of  the  curve  of  the  objective 
function. 
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3.3  A  Chemical  Manufacturing  Problem  with  Advertisement 

The  >'.o..;ol 

Figure  34  represents  a  chemical  manufacturing  process  and  stages  1 
and  2  represent  two  reactors.   The  raw  material  entering  the  first  reactor 
is  a  mixture  of  A  and  B.   After  the  second  stage,  the  product  A  and  product 
B  are  separated,  as  is  the  remaining  raw  material,  product  C.   Product  B 
is  the  more  valuable  of  the  three  products  and,  to  enhance  its  sale,  it  ' 
has  to  be  advertised.   Also,  to  meet  the  fluctuations  in  its  demand,  a 
certain  amount  of  inventory  has  to  be  kept.   It  shall  be  assumed  that  the 
demands  for  products  A  and  C  are  unlimited. 


A(t) 


SALES 


V   A,C 


Fig.  3A 
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Let  x~,  and  y~  represent  the  concentration  of  A  and  B  in  the  original 
raw  material  before  it  enters  the  first  stage  or  reactor.   Similarly, 
let  x  ,  y   and  x  ,  y   represent  the  concentrations  of  A  and  B  before  and 
after  the  second  stage,  respectively.   To  bring  about  this  reaction, 
temperatures  T..  and  T_  have  to  be  applied  to  the  two  reactors.   The  re- 
actions in  the  reactor  can  be  represented  by  the  following  equations: 
Let  q  =  flow  rate 

v  =  volume  of  the  first  reactor 

v~  =  volume  of  the  second  reactor. 

Then, 

dxl 
Vl  dT~  =  q(x0  '  Xl)  "  vi  Kai  xi  (79) 

dyl 

vi  dT"  =  q(yo  "  yi}  "  vi  Kbi  yi  +  vi  Kai  xi  C80) 

dx2 
V2  dT  =   q(xl  "  X2}    "  V2  Ka2  X2  (81) 

dy2 
V2  dT  =  q(yl  "  y2}    "  V2  Kb2  y2  +  V2  Ka2  X2  (82) 


wnere 


Kan    =  G     exp    (-Ea/RTj 
i  a  1 


Ka2  =  G  exp  (-Ea/RT  ) 
Kb1  =  G  exp  (-Eb/RT  ) 
Kb2  =  Gb  exp  (-Eb/RT2). 
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This  completes  the  production  part  of  the  system.   Now  consider  the 
;;tory.   The  rate  of  change  of  inventory  is  the  difference  between 
rate  of  production  of  B  and  its  rate  of  sale.   If  I(t)  represents  the 
inventory  at  time  t,  then 


dl(t)         „  „,  N 
JT2-   =  ^2  "  Ca  K(t) 


(83) 


The  sales  equation  is  assumed  similar  to  the  problem  in  Para.  3.2. 


dK(t) 
dt 


=  [C  +  A(t)]  •  K(t) 


1  - 


K(t) 

N 


(84) 


Equations  (79)  through  (84)  represents  the  performance  equations  of  the 
whole  system  under  consideration. 

This  problem  has  six  state  variables,  namely  x, ,  y, ,  x  ,  y„,  I(t). 
K(t)  and  three  control  variables  namely  T.  ,  T~  and  A(t) . 

The  profit  function  can  be  formulated  as: 
Profit  =  (sales  revenue  from  A,B,C)-(cost  of  holding  the  inventory  for 
B)  -  (cost  of  advertising  for  B)  -  (cost  of  production) 
Sales  revenue  from  A,  B  and  C  is  =  C,C  K(t)  +  C?qx  +  CLq  (1  -  x   -  y  ) 

where,  C.  ,  C_,  C,  represent  the  unit  sales  prices  for  A,  B,  C  respectively, 

2 
Cost  of  holding  the  inventory  of  B  =  CT  (Iw  -  I(t))   where  I__  is  the 

I   M  M 

capacity  of  the  warehouse  and  C  =  inventory  carrying  cost. 

2     2 
Cost  of  advertising  =  C  A  (t)  K  (t) . 

A 

Cost  of  production  comes  from  the  fact  that  the  two  reactors  have  to 
be  supplied  with  heat  energy  in  order  to  obtain  the  desired  temperature. 
Let  C   represent  the  cost  of  raising  the  reactor  temperature  by  a  unit 
degree.   Then  the  cost  of  production  becomes 


"  CT{(Tlm  "  ^  +  (Tl-T2)2} 


where  T..   is  the  temperature  of  the  entering  raw  material.   Thus  the 
function  to  be  maximized  is 
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J  =  /   C..C  K(t)  +  C0qx. 

U        1  q  2  2 


C3q(l-x2-y2) 


CT[I  -I(t)] 
I  m 


-  CAA2(t)K2(t)  -  CyCT^-T^2  +  (T^T^2]}- 


£  -  ClCqK(t)  +  C2qx2  +  C3q(l-x2-y2)  -  C^-lCt)]2 


-  CAA2(t)K2(t)  -  QjKT^)2  -  (TrT2)2]   . 


(85) 


Recursive  Relations 

The  necessary  relations  for  the  second  variation  can  be  obtained 
in  the  following  manner.   The  various  derivatives  can  be  obtained  as 
follows 


3x 


3x. 


3y- 


3J 

3x, 


3J 

3y2 

3J 
31 

3J 

[  3K   , 


q(C2-C3) 


"C3q 


2CT(I  -I(t)) 
1  m 


C1Cg-CAA  (t).2.K(t) 


q  a 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3x 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-2CI 

0 

0 

0 

0 

0 

0 

-2CAA2(t) 

a    J 
363x 


0  0 


0  0 


_3J 

36 


2CT(Tlm-T1)   +   IC^I^T^ 


2CT(T1-T2) 


-2CAA(t)K    (t) 


36' 


-2C„ 


2C, 


-2C. 


-2C.K    (t) 
A 
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3ff 

3x 


BM      Ga-EAT 


BM, 


(q/V, 


BM, 


q/V2 

Ga-EAT, 

BM, 


-C       BMC 

q  5 


where, 


BM1  =  -  (q/Vx)  -  Ga-EAT.^ 


EAT1  ■  «P<"  W? 


BM2  =  -  (q/V1)  -  Gb-EBT  , 


EAT2  "  «*<"  W? 


BM3  =  -  (q/V2)  -  Ga-EAT2,  . 


EBT1  "  ***<-  W? 


BM4  =  -  (q/V  )  -  Gb-EBT^ 


EB 
EBT2  -  «p(-  — ) 


BM5  =  (C+A(t))  -  Gb-EBT2 
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L3x    J 


BM, 


Ga-EAT  BM  0  0  0 


q/V, 


0  BM, 


0  q/V2     Ga-EAT2     BM^  0 


-C 


BM 


5  J 


af ' 

36 


FT. 


FT  0  0  0  0 


FT-        FT^  0  0 


0  FT 


5  J 


vhere 


FT,    =  - 


Ea      „  RT1 

~~T  GaXle 
RT^ 


Eb 


FT     = 
2 


Eb 


RT 


GJ^ 


RT 


2     b*7! 


1  Ea      _ 

+  5* G^e 


Ea 
RT, 


Ea 


Ea 


FT,  =  -  — -r  G  x  e     RT? 
3  RT^     a2  2 


105 


FT^  = 


Eb 
RT 


Eb 
RT. 


2  °by2e 


Ea 
RT. 


FT5  =  K(t)  [1  -  ^|i] 


82f. 


3x 


3"f, 


8x 


32f. 


8x 


2      2 

til    lis. 


3x 


3x 


0 

0 

0 

0 

0 

'0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

'  0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

32f 
3  f6 

0 

0 

0 

0 

0 

0 

3x2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-  |  [C+A(t)] 
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3x 


r  BM1 

0 

0 

0 

0 

.      Oj 

32£. 


3x  80 


0 

0 


3x 


r 

3a. EAT. 

BM2 

0 
0 
0 
0 


32f. 


8x  de 


-G 


0     0     0     0 

0     0     0     0 
0     0     0     0 


3f. 


3x 


q/V, 


-q/V2-Ga.EAT2 


32f. 


ix  3( 
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3f 


0 
q/V. 


Ga-EAT, 


BM, 


9x  36 


o 

-G. 


af, 

3x" 


0 
0 
0 

q 
o 

-c 


32f. 


3x    86 


'   0 

0 

0 

0 

0 

0   ' 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3x 


0 

0 
0 
0 
0 


[i  -^-Hc-fA(t)] 


2r 

r  o 

0 

0 

0 

0 

0 

f. 

6 

3x   36 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

[1 

2K(t) 

N 


where, 
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°l" 


EA     -  .     EA« 

— -  Ga   exp(-  — ) 

RT?  OT1 


EA 

RT' 


G_   =   -  — g  Ga    exP("  — 2 


EA 


EB      _  ,      EB    > 

G2  =  -  —  Gb  exp(-  g-) 


EE      „,  /      Eb    n 


5 

36 


ha 


Ea 
RT, 


G  x.e 


RT 


2      a  1 


32f. 


86 


GlXllT  2; 

X  X  Al        RT^ 


3f, 

36" 


Eb 

- 

Eb 

RT^ 

GbyiG 

RT 

+ 

Ea 

2 
RT^ 

G  x.  e 

a  x 

Ea 
RT 

32f, 


36 


/2  Eb     s 

"  G2yl(T"  "  "~2   } 
^   x   xl        RT^ 


/      2  Ea     v 

+   G  x   (-  —  ♦  —  ) 

X   -1        Xl        RT^ 
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3f 
38 


3  _ 


Ea 


Ea 
RT, 


G  x^e 


RT 


2     a  2 


32f 


3  _ 


30 


„    I      \/2          Ea   y 
G   (x    )(- -) 

J  2       RT* 


36 


Eb 
RT 


Eb 
RT, 


>   V2e 


Ea 
RT, 


a2f, 


36 


/2  Eb    x 

GUy2(~  '  ~) 
d  T-        RT* 


i      2  Ea    \ 

G3X2(-T~+~T) 
J  2       RT* 


3f 
3? 


>   _ 


0 
v.       j 


32f, 


36 
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Hi 

30 


K(t)(l  -  ^-) 


32f, 


30 


0     0 
0     0 


P  = 


11 


21 


31 


Ul 


51 


61 


21 


22 


32 


U2 


52 


62 


31 


32 


33 


U3 


53 


63 


Ul 


U2 


U3 


UU 


5^ 


6U 


51 


52 


53 


5U 


55 


65 


61 


62 


63 


6k 


65 


66 


Now  the  expressions  for  R_,  S_,  T  which  are  required  for  obtaining  the 
second  variational  Equations  (19),  (28)  and  (U8)  may  be  determined. 
Thus 


p  ..   52J    3f ' 
-   39  3x 


32f. 


il_p  +  y  z  — i 

36  -   A  Zi  30  3 


i=l 


3x 


Let 


R  = 


10 


11 


13 


llU 


ll6 


IT 


12 


15 


V18 


Ill 


where 


Rl  =  Pll  m  +  P21  ^  +  Zl°l  "  Z2°l 


R  =  P   FT3  +  Pj.  FTU 


R3  =  P6l  FT5 


R^  =  P21  FT1  +  P22  FT2  +  ZgG2 


R.  =  P32  FT3  +  P^g  FTl4 


R6  =  Pg2  FT5 


R_  =  P31  FT1  +  P   FT2 

Rg  =  P33  FT3  +  PU3  FTl*  +  Z3G3  -  ZUG3 

R9  =  P63  FT5 


R   =  PUl  PT1  ♦  P^  FT2 


Rll  =  Pl.3  FT3  +  Tkk   FTU  +  \°U 


R12  ■  p6i,  ra5 


R13  =  P.   FT1  +  P   FT2 


RlU  =  P53  TO  +  P51.  m 
R15  =  P65  PI5 


R16  =  P6l  FT1  +  p62  ™ 


R17  =  P63  m  +  P61.  m 


From  Equation  (U6) 

s  "  ae"  +  *  zi  aF" 
—  1=1    — 
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S  = 


S 


2ct(ti»-ti)  +  **<*i-*a> ♦  ziGixi +  (Vi 

2CT(Tl-T2)  ♦  z3x2G3  ♦  zu(y2Gu  -  x^) 
z6  K(t)    [1  -  ^-]    -   2CAA(t)   K2(t) 


GlXl'    *2 


From  Equation   hj 


T  =  i-^L+     y     z    L 

-  2         *•       i         2 

3G  i=l  30 


Let 


T  = 


(DTI         DTU  DT?' 


DT2  DT5  DT9 


wDT3  DT6  DT9 


wnere 


DTI  =  -   z  G  x      (—  -  -^r)   -   G  z  y      (— 
111       Tx        RT2  2   2   1      Al 


Eb    x 
2 


f        2        Ea 
22G1X1    (_T7+^2 


DT2   =   -  2CP 


DT3  =   0 
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DTU  =  2CT 


„„^      _„      _     /2    Eel  \      _     /2    Ed  \ 
DT5  -  -  89,  -  ,3G3x2  (5-  -  -£>  -  z^y2  (5-  -  -|) 


/    2    ,   2    Ea  v 

-  G3v2  (--+  (-~  +  — ) 

J        2       2   RT^ 


DT6  =  0 
DTT  =  0 

DT8  =  0 


DT9  =  -  2CA  K2(t) 
dz 


The  equations   for  —  ,   i.e.  Equation   (19),   take  the  following  form: 

CL  \f 

dz 

— -  =  -   z.BMl  -   z_G     EAT1  -    (q/V.)    z,  (86) 

dt  l  d  a  ^       o 


dz0 

^=  -   z2BM2   -   zu(q/V2)  (87) 


dz 

^2.  =  -  q(C2-C3)    -   z  BM3  -   G&  EAT2   z^  (88) 


dzi 

-rr=-  =  C_q  -  BMl+   z,    -   zc   q  (89) 

dt  3  ^5 


dz 

^=2Cl    (IM-Kt))  (90) 


az6  ? 

~~-  =  -   C,C     +   C     A   (t)    2K(t)   +   zcC     -   z,BM5  (91) 

dt  1   q  a  5   q  o 


Let 
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R'TR 


RTR1 

RTR2 

RTR3 

RTR4 

RTR5 

RTR6 

RTR2 

RTR7 

RTR8 

RTR9 

RTRIO 

RTR11 

RTR3 

RTR8 

RTR12 

RTR13 

RTR14 

RTR15 

RTR4 

RTR9 

RTR13 

RTR16 

RTR17 

RTR18 

RTR5 

RTR10 

RTR14 

RTR17 

TRT19 

RTR18 

RTR6 

RTR11 

RTR15 

RTR18 

RTR20 

RTR21 

Equation  (28)  now  becomes 


dP 


n 


- =  -'2Pin    BM1  -   2Pon    Ga  EAT1  -   2PQ1    (q/V0)   +  RTR1 

at         11         21  61  2. 


(92) 


dP 


21 


=  -  P21  (BM1  +  BM2)  -  P22  Ga  EAT1  -  (P32  +  P^) (q/V2)  +  RTR2 


(93) 


dP 


31 


-r-^  -  -  P__  (BM1  +  BM3)  -  P..  Ga  EAT2  -  P._  Ga-EATl  -  P_„  (q/V0) 
at        Jl  41  JZ  jj     / 

+  RTR3 


(94) 


dP 


41 


-  -  P..  (BM1  +  BM4)  -  P._  Ga  EAT1  -  P.,  (q/V0)  -  P,,  q  +  RTR4 
at         41  4Z  4J      Z      jl 


(95) 


dP 


51 


=  -  P_.  BM1  -  Pco  Ga  EAT1  -  Pco  (q/V„)  +  RTR5 

at      ji      dz  _>j    z 


(96) 


dPfil 

-j-^  =  Pcn  C   -  P,.  (BM1  +  BM5)  -  P,0  Ga  EAT1  -  P,«  (q/Vj  +  RTR6 
dt      51  q    61  62  63   n   2 


(97) 


dP 


22 


■j^-   =  -  2P22BM2  -  2PA2  (q/V2)  +  RTR7 


(98) 


dP 


32 


=  -  P32(BM2  +  BM3)  -  PA2  Ga  EAT2  -  P ^    (q/V,,)  +  RTR8 


(99) 


+  RTR21 
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dP 


U2 


=  -  PU2  (BM2  +  BMli)  -  Puu(q/V2)  -  P^  q  +  RTR9 


(100) 


dPS? 

-^L   =  -  P   BM2  -  F^    (q/V2)  +  RTR10 


(101) 


dPfip 

at   =  P52  CQ  "  P62  (M2  +  BM5)  "  ?6k    (q/V2)  +  RTR11 


(102) 


dP 


33 


=  -  2P   BM3  -  2P,   Ga  EAT2  +  RTR12 


(103) 


dP, 


43 


dt 


=  -  P,   (BM3  +  BMU)  -  Pn  Ga  EAT2  -  P„»q  +  RTR13 


U3 


kk 


53 


(10U) 


dP 


53 


dt 


=  -  P        BM3  -  P5l+   Ga  EAT2  +  RTRlU 


(105) 


dP 


63 


=  Pco  Cq  -  P,_  (BM3  +  BM5)  -  P^,  Ga  EAT2  +  RTR15 
dt      53       63  oh 


(106) 


dP 


kk 


dt 


=  -  2P^  BMU  -  2P  ^   q  +  RTR16 


(107) 


dP  , 

—2-*-  =  -  P  .  BMU  -  P„  Q  +  RTR1T 
dt        5k  55 


(108) 


dP,-. 

— —  =  Pci  Cq  -  P.,  (BMU  +  BM5)  -  P,c  q  +  RTR18 
at      5^       64  65 


(109) 


dP 

-r-^-  =  2C_  +  RTR19 
at       I 


(110) 


dP^, 

=  Pcc  Cq  -  Pr_  BM5  +  RTR20 
dt      55   *    65 


61 


(111) 


dP.^ 

^^  =  2z6  [C  +  A(t)  ]  +  2Ca  U   {t))   +  2P^  Cq  -  2P66  BM5 


(112) 
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Equation  (^8)  is  given  by 


dt 


=  R'  T"1  S  +  R'  T"1  (~ -)  &  -  (j|4  £ 


where  ^  is  six  dimensional.   Here  all  the  terms  were  obtained  by  the 
matrix  multiplication.  The  new  control  is  calculated  as  given  by 
Equation  (50),  i.e. 


'    -'_         | 

T2 

.  A(t)   . 

(J+l) 


x2 

I  A(t), 


(J) 


-      ^-36 


c  T  "(S  +  ■■£-  &) 


(J) 


-  (T^RjJ  •  (x(J+1-x(^) 


Table  18 


Numerical  Values  of  the  Constants 
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Set  h 

1  1 

q 

-  60 

c 
q 

=  1.0 

Im 

=  20.0 

C2 

=  0.0 

CA 

=  0.0002 

EA 

=  18000.0 

yi 

=  0.430 

Set  // 

1  2 

q 

=  60. 

c 

q 

=  1.0 

Im 

=  10. 

C2 

=  0. 

CA 

=  0.01 

EA 

-  18000 

vl 

= 

12.0 

c 

» 

1.0 

lm 

" 

340. 

C3 

= 

0.0 

CT 

= 

0.0005 

EB 

= 

30000.0 

Ga 

- 

0.535x10 

11 


=  12.0 


yx  =  0.43 


vl 

J-  *-  •  \J 

c 

= 

1.0 

lm 

■ 

340. 

C3 

= 

0 

CT 

= 

0.001 

EB 

= 

30000 

Ga 

zr 

0.535x 

10 

v2 

= 

12.0 

N 

= 

100.0 

Cl 

= 

5.0 

CI 

a 

1.0 

R 

= 

2.000 

XI 

= 

0.530 

Ga 

- 

0.461x10 

18 


11 


v2 

N  = 

ci  = 
ci  = 

R  = 

XI  = 
Ga  = 


12.0 

100 
5.0 
1.0 
2.0 

0.53 

0.461x10 


18 


Set  //  3 

Same  as  Set  //  2 
except     Im  =  20. 


'  CT  =  0.0005 


and        K(to)  =  lt0 
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Discussion 

This  particular  problem  reveals  how  the  theoretical  attractiveness 
the  second  variation  method  is  more  than  offset  by  both  the  complexity 
and  by  the  number  of  equations  to  be  integrated. 

In  this  problem  (6+1)  or  seven  equations  are  to  be  integrated  in  the 
forward  direction  and  (6+6(7) /2+6)  or  33  equations  in  the  backward 
direction.   In  addition,  the  calculations  of  R,  S,  and  T_  are  in  terms  of 
matrix  multiplications  and  T   has  to  be  calculated  at  each  step  of  the 
integration  in  Equations  (28)  and  (48) . 

This  program  was  tried  with  three  different  sets  of  numerical  values 
which  are  shown  in  Table  18.   These  values  were  taken  from  the  solution 
of  the  same  problem  first  by  variation  and  quasilinearization  respectively, 

This  problem  was  found  to  be  unstable  as  far  as  its  solution  by 
the  second  variation  is  concerned.   With  all  the  various  values  tried, 
the  program  could  make  a  complete  iteration.   However,  it  fails  in  the 
backward  integration  of  the  second  variational  equations  because  of 
exponential  overflow. 
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U.      CONCLUSION 

The  second  variation  method  has  been  shown  to  be  a  fairly  useful 
tool  for  attacking  the  complex  practical  optimization  problems  involving 
a  fairly  large  number  of  variables.  The  convergence  is  very  fast,  pro- 
vided the  initial  or  starting  guess  is  sufficiently  close  to  the  optimal 
trajectory.   This,  however,  becomes  more  and  more  difficult  when  more 
than  one  control  variable  are  involved.   In  that  case,  the  number  of  combin- 
ations that  could  be  used  as  the  starting  trajectory  is  quite  large  and 
makes  the  initial  guess  a  difficult  task.   This  can  be  overcome  by  using 
the  first  variation  method  for  the  first  few  iterations  and  then  switching 
to  the  second  variation  method.   This  combination  provides  rapid  conver- 
gence to  the  optimum  from  almost  any  realistic  starting  trajectory.  The 
theoretical  attractiveness  of  this  method  is  removed  by  its  disadvantages 
like  the  guess  of  the  initial  trajectory  for  the  state  variables  in  ad- 
dition to  that  of  control  variables.   Also  the  number  of  equations  and 
their  complexity  make  the  use  of  this  technique  tedious. 

The  first  variation  method,  of  which  the  second  variation  method 
is  a  natural  evolution,  should  be  used  in  combination  with  the  second  vari- 
ation.  The  first  variation  method,  unlike  the  second  variation,  will 
approach  optimum  from  almost  any  realistic  starting  trajectory.   The  re- 
sults of  the  first  variation  method  could  then  be  used  as  the  starting 
trajectories  for  the  second  variation.   In  this  way,  the  convergence  problem 
of  the  second  variation  can  be  partly  overcome.   This  combination  provides 
a  rapid  convergence  from  almost  any  realistic  starting  trajectory  for 
most  engineering  problems.   While  evaluating  the  merits  and  demerits  of 
this  technique,  it  should  be  borne  in  mind  that  no  single  optimization 
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technique  is  suitable  for  all  classes  of  problems  that  will  be  encountered. 
Each  technique  will  be  most  efficient  only  for  a  particular  type  or  types 
of  problems.   It  is  left  to  the  decision  of  the  engineer  to  select  any 
one  or  a  combination  of  these  techniques  for  the  problem  he  is  facing. 
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7.   APPENDIX 

7.1  Computer  Program  for  the  Inventory  Model 

7.2  Computer  Program  for  the  Inventory  and  Advertising  Model 

7.3  Computer  Program  for  the  Chemical  Manufacturing  Problem  with 
Advertisement 


s 

1C 


$JI 
c 


=CHECK,TIME=9,PAGES=200,L      >0    126 

[ATI   .  -i  ,    ADIENT  rECHNIQUI 

Pi  KENS  I  ON  Pl(  1  50)  ,P2(  150)  ,P3(15u)  ,P4(  150),  UX1  (  15Q),X1  (  150)  ,(.■!(  15  0) 
Olf  "   I  N  Q2(  150)  ,U  3(  150)  ,UM  150)  ,0X2(  150)  f  X2(  150) 

SION  Z(150),M(15<  ),T(150) 
DIH1  \510N  RSM  150)  ,KT  (  150)  ,K(150) 


100 


(HA 


1 


I   1  L( 
T  (1H-, 

P 


150)  ,  U(  lr><  )  ,X(  15  0)  ,P(  150)  ,TX(  150)  ,Q(  150)  ,  S(  150) 
•NO.  XI  I 

0  TI'.ETA 


I   '  I 


10  1  F0RMA1  (IH-.tlM  ,13,2X,E15.8,2X,E15.8,2X,E15.8,2X,E15.8,2X,E15.8,2 

IX, I  15.8, IX, E 15.8) ) 
10/  FORHAT  18F9.7) 


103 

104 


FC 
I  M).1 


1   (  I  H  ,  9  (  F  7  .  1)  ) 

T  ( I H-, • **********************♦*♦*********#******    tTERATI 


11 
12 
13 


I TMAX=58 
.1 

103, P 


A,CP,XM,CI,D,fJ,A,XHl),X2(l) 
PA,CP,XMfCl,0fB?AfXl(l),X2(l) 

101 


14 
15 
It 


PRINT 
CP=.C 
DO  12 


103, 

01 
1  =  1, 


17 
it 
IS 


12 


n  1  )  = 

KC I  >- 

DO  I 


). 

5. 

N=l,  ITM<\X 


2C 

2  1 
22 


LCN>  = 

DO  6 
"'([)  = 


N 

[  =  1,1 

I 


01 


2  3  Pl( I  )=D*( T (  I  )-A- 

24  P2(  I  )=D*(T(  I  )-A- 


B*D*( 1-1 ) ) 

B*(D*( [-!)«■. 5*D) ) 


FOR**RD  INTEGRATIlj   II   /I 


25 
26 
27 


M(  I  ) 
P4(  [  ) 
DXK  I 


=  P2(  I 

=D*(  r 

)=(  1 . 


) 

( I )-A- 
/5. )*( 


B*(D*( 1-1 >+D) ) 

Pl(  I )+2.*P2( 1 )+2.*P3(  I  )+n4(  I  )  ) 


2£  Xl(  l  +  l)  =  Xl(  I  )+DX 

25  )1(1)=D*(CI*(XM- 


33 
34 

35 


111) 

**$***************    FORWARD  INTEGRATION  OF  X2 

XI  (  I  )  )*(XM-X1  (  I  )  )+CP*EXP(  (PA-T (  I  )  )*(PA-T  (  I  )  )  )  ) 


3C  g2(l)=D*(CI*(XM-Xl(I)+.5*Pi(I))**2+CP*EXP((PA-T(I))**2)) 

31  i (  I  )=D*(C[*(XM-XK  1  )+.5*P2(  T  )  )**2*CP*EXP(  (PA-T(  I  )  )**2) ) 

32  04(l)=O*(C[*(XM-Xl(lH-P3(I))**2  +  CP*EXP((PA-ni))**2)) 


HX2(  I  )=(  1./6.  )*(Q1(  I  )+2.*Q2(  I  )+2.*Q3(  I  )+Q4(I  )  ) 

X2(  1  +  1  )=X2(  I  )+DX2(  I  ) 

TX(  I  )=2.*CP*(EXP(  (PA-T(  [  )  )**2)  )*( i.*2.*(  (PA-T (  I )  )**2)  ) 


36 
37 

3e 


S(  I  )  = 

00  / 
K=IC1 


2.*CP*(PA-rm)*(Exp(iPA-rm)**2)) 

1=1,100 


35 

4C 
41 


ZIK1 

PI  I'M 

0(1.1 


)  =  0. 
)=0. 


42 


l  =  - 


.♦C1MXM-XKK-H-1  )  ) 


BACKWARD  INTEGRAriO\  OF  / 


43 
44 
45 


ZM2=Z 

*=Z 

ZM4=Z 


Ml 

Ml 
Ml 
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46  »    Z(K-I)=Z (K+l-I )+(l./6.)* 

Casfc**^*:**^^*  ***************** 


(ZM1+2.*ZM2+2.*ZM3+ZM41 

*******BACKWARD    INTEGRATION    OF    DP/D1 

***  DP/Cr=-2Cl+ (P) SQUARE ( IX  ) 


47 

48 
4S 


PM1=-D*(-2.*CI+TX(K+1-I )*(P( 
PH2=-D*(-2.*CI+TX(K+1-I)*( (D 

PM3=-D* ( -2. *CI+rx(K +!-!)*( (P 


K+l-I )**2) ) 

(K+l-I J+PM1/2. >**2I ) 

(K+l-I  )+PM2/2.  )**2)  ) 


5C  PM4=-D*(-2.*CI+TX(K+l-I )*( (P 

5  1  P(K-I)=P(K+1-I)+<1./6.)*(PM1 


(K+l-I  )+PM3)**2)  ) 

+2.*PM2+2.*PM3+PM4) 

******BACK^ARD     INT  EGRAT  ION 


Q$$$$*:«::j£4:*$$*£#**$**************** 

52  l=-D*(  (P(K+1-I  )*(  Z:  (  K+l-I  )  + 

53  QM2=-D*( ((P(K+l-I)+.5*PMl)*( 


******  DQ/DT=(  PIZ  +  O-S)  )/TX 

0(K+1-I)-S(K+1-I))  )/TX(K+l-I  )) 
Z(K+1-I)+.5*ZM1+Q(K+1-I)+.5*QM1-S(K+1 


54 


II  )  )  )/ 

QM3*- 

II  )  )  )/ 


T  X  (  K  + 
D*<  (  ( 
TX(K  + 


1-1  )  ) 
PIK+l-I 
1-1  )  ) 


)  +  .5*Pf-'2  )*(Z  (K+l-I  )  +  .5*ZM2  +  Q(  K+l-I  )  +  .',  -S(K+1 


55 


56 


0M4=-D*( UP(K+l-I)  +  PM3)*(Z(K+l-I)+ZM*  +  Q(K  +  l-I ) +QM3-S (K+l-I ) ) )/ 
l  +  l-I  )  ) 
7    Q(K-I  )  =  0( K+l-I  )  +  (l./6.  )*(0M1+2.*QM2+2.*QM3+GM4) 


57 
58 
5S 


DO    8    1=1,101 

ri(i)  =  T(i)-EP*((7(i)-S{i)+g(i))/rx(i))-(P(i)/Tx(i))*(xi(i)-)((i)) 

P'UNT     104,  LIN) 


6C 
61 
62 


PKINT 
PKINT 
DO  -7 


100 
101, 
1  =  1,1 


(M(I),Xl(l),X2(I),?(I),P(I),G(I),TKI),rx(I),I=l,101) 

01 


63 

64 
65 


T(  I  J  = 
9  X  (  I  )  = 
I  CHNTI 


Tl(  I  ) 
XI  (  T  ) 

NUE 


66 
67 


STOP 

END 


4 

c 

6 
7 
6 

1C 

11 
12 
12 
14 
15 

16 


17 
1F 

i<; 

2C 
21 
22 
23 
24 
25 

26 

27 

26 

2S 

3C 

31 

32 
2  2 

34 
35 
36 
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FXIFRNAL  FCT,OUTP  128 

DIMENSION  PRMT(iO),CERY(lO) , AUX(8,  12  )  , YO ( 2 ,  10  1  )  , Y  1  ( 9 ,  10  1  ) , Y (  10 ) , AT 
KlOU  ,ATNFW(  101) 
COMMCN  Yl,/\T,ATNEhtYnfAlB,C,AN,FlCI|P[,CA,KfNK,FP,RlfR2ib,T,J,* 

C   THIS  PROBLEM  HAS  2  STATE  VARIABLES  AND  1  CONTROL  VARIABLE 

C  —  THL  STATE  VARIABLES  ARE  THE  INVENTORY  AND  THE  SALFS 

C  --  THE  CCNTROL  VARIABLE  IS  THE  ADVERTISEMNT  MACE 

C 

C   — 
ICO 

101 

102 
103 

104 

105 


C   VA 
2C0 


201 


2C3 


2C2 


FORPA 

REAC 

FO.tMA 

REAC 

FORMA 

FORM 

l=',FE 
FORMA 

1=* ,F8 
FORMA 
PRINT 
PUNT 
PRINT 
PRINT 

LUES  C 
FORMA 

101  GR 
PRINT 
FORMA 
REAC 
REAC 
REAC 
FORMA 
PKINT 
FORMA 
PRINT 


T  (9F 

100, A 

T  (21 

101, N 

T  C 

T  CO 

.3) 

T  (  '0 

.3) 

T  (  «0 

102 

103, 

104, 

105, 

F  THF 

T  CO 

ID  PO 

200 
T  (12 
201,  ( 
201,  ( 
201,  ( 
T  (  '0 

203 
T  (  IH 
202, 


RCA0  THE  VARIOUS  CONSTANTS   

8.3) 

,B,C,AN,F,CI,PI,CA,EP 
4) 

Kf ITMAX 

FOLLOWING  VALUES  OF  THE  VARIOUS  CONSTANTS  ARE  READ  tN« ) 

B=« ,F8. J, • 


F=»,F8.  U 


CI=' ,F8.3, • 
NK=«  ,  14,  ' 


C=» ,F8. 3, 


PI=f ,F8. ), ' 


AN 
CA 


ITMAX=« ,14) 


A,B,C,AN 
F,CI,PI,CA 

EP,NK, ITMAX  

STATE  ANC  CONTROL  VARIABLES  AT  101  GRID  POINTS  

THE  FOLLOWING  VALUES  OF  STATE  AND  CONTROL  VARIABLES  AT  1 
INTS  ARE  READ  IN'  ) 


F6.3) 

Y0( 1,NM) ,NM=1,NK) 
YC(2,NM) ,NM=1,NK) 
AT(NM) ,NM=1,NK) 

NO.  INVENTORY 


SALES 


ADVT. • ) 


C   VA 


DO  2 
FORMA 

I  NO. 
PRINT 

RICUS 
PRMT  ( 
PRMT( 
PRMT( 
PRMT( 
NUlf  = 
DERY  ( 
DERY  ( 
00  4 


-,  ( IH  ,  I3,4X,3(F6.3,5X)  )/) 
(NMf Y0( 1,NM),Y0(2,NM) ,AT(NM) ,NM=1,NK) 

^AIN  cn  LOOP  FOR  ITERATIONS    

IJ=l,  ITMAX 

j     (iQ*********************************************   ITERATION 
^14,*  *************************************  •) 
J,  IJ 

PARAMETERS  FOR  FORWARC  INTEGRATION   

l)=0. 

2)  =  l. 

3)  =  .0 
4)=.0 
3 

1)=ND 
11*1. 
1=2, N 


1 
0  1 

IM 

/DFRY( 1 ) 
IDIM 


37 

3e 

3S 


4C 
41 

42 

A3 


44 
46 


41 
4€ 
4<5 


5C 
51 
52 


53 
54 


55 
56 
57 


56 
5S 
6C 


61 
62 


DERY(  I  )=DERY(  1) 

Y  (  I  )  =  .  ? 

Y(2)=.2 

Y( 3)=C. 

K=i 
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CALL  RKGS  FOR  THE  FORWARD  INTEGRATION 

CALL  RKGS(PRPT,Y,DERY,NDIM, IHLF , FCT , CU TP , AUX ) 
VARIOUS  PARAMETERS  FOR  THE  PACKWARC  INTEGRATION 

PRMTd  )  =  1. 

PRMT(2)=0, 

PRMT(3)=-.01 

PRMT(4)=.01  

NDIM=7 

DERYU )=NDIM 

DERY(  l)  =  l./DERY( 1 )  

DO  5  I=2,NDiy 

5  DERYI  I)  =  DERY(  I) 

DO  6  1=1,7  

6  Y( I )=0. 
K  =  2 


—  CALL  RKGS  FOR  BACKWARD  INTEGRATION 


CALL  RKGS(PRMT,Y, DERY,NDIM,  IHLF f FCT ,CU TP , AUX ) 
DO  7  L-lfNK 

YO(  1,L)=Y1(  1,L)  

YO(2,L)  =  YH2,L) 
AT(L)=ATNEW(L) 

CONTINUE  

STOP 
END 


6<; 

7C 
71 
12 


72 
74 
11 
It 

77 
76 
7S 

ec 

81 

ei 

82 


84 
85 

ee 

67 


c 
c 
c 


c 
c 

c 


SUt  RC 

DIMEN 

l(l<  1  ) 

CCM^C 
CEPFNI3IN 
CF  THIS 
INTEGRA! 

IF  (K 
—  THIS 

IF  (  X 

J=L 

Y( 1 } ,Y(2) 

INUCPt  M 

PROGRAM 

il  DERYl 

OERYI 

DE 

RETLR 


SE 

10  IF(X. 

N=NK 
1?  Rl  =  Y< 

R2=-2 
IN) 

S=-2. 

T  =  -2. 

DERYI 

DERYI 


UTINIE  FCTIX,  Y,  CFKY)  130 

SIQN  PRMT I 10),DERY( 10), AUX( 8,12),Yn(2,lCi»,Yl(9,lOl),YlU  ),  J 

,ATNEW( 10L) 

N  Yl,AT,ATNE*,YC,AfB,C,AN,FtCI,PI,CA,K,NK,FPfRl,R2,S,T,J,N 

G  ON  THE  VALUE  OF  K,CITHFR  THE  FIRST  PART  OR  THE  bECOND  PART 

'OUTINE  IS  USED  FOR  THE  FOKWARD  AND  THE  BACKW/. 
ION  RESPECTIVELY. 
.EQ.2)  GO  TO  1C 

PART  OF  SUBROUTINE  IS  FOR  FORWARD  INTEGRATION  ONLY   

.NE.O)  GO  TO  11 


DENOTE  X  ANC  G  IN  THE  PROBLL."-   

ENT  VARIABLE  T  IN  THE  ORIGINAL  CCNS.  IS  UENCTED  BY  X  IN  THE 

1  )  =  A*B*X-Y( 2 ) 


2)=Y(2)*(C+AT(J))*(1.-Y(2)/AN) 

H«YI 2)*F-CI*( (PI-YI  I)  )**2)-CA*YI2)*( ATIJ)**2) 

N 

FIRST  PART  FOR  FORWARD  INTEGRATION  ENDS   

COND  PART  FOR  EACKwARD  INTEGRATION 

NE.l)  GO  TO  12 

4)*YK2,N)*(  l.-YK2,N)/AN)    

.*CA*AT(N)+Y(5)*Yl(2,N)*(l.-Yl(2fN)/ANUY(2)*(l.-2.*YH2,M/A 

*CA*Y1(2,N)*AI (N)+Y(2)*Y1(2,N)*{1.-Y1(2,N)/AN) 


C 

c 
c 
c 


DERYl 

UFRYI 

DERYI 

1 (2,N  ) 

DAC 

TF 

CF  IN  TF 

THIS  PR 

DERYI 

DERY  I 

l)*l  1. 

REILM 


*CA*Y1(2,N) 

1  )=-2.*CI*lPI-Yl I 1,N)  ) 

2)=-F  +  CA*(AT(N)**2)+Y( 1 )-Y I  2  )  *  I C  +  AT  I N  )  ) *  I  1 .-2 . *Y 1  I  2 , N ) / AN ) 

BACKWARD  INTEGRATION  OF  DP/DT  3  EQUATIONS   

3)=2.*CI+(R1**2)*T 
4)=Y(3)-Y14)*IC  +  AT(N))*U.-2.*Y1  I  2  ,  \  )  /  AN  )  +R  1  *K2*  T 
5)=(2.*Y(2)/AN)*(C+AT(N))*2.*YU)-2.*Y(5)*(C+AT(N»»*(l.-2.*Yl  1 
/AN)+(R2**2 )*T 

KWARD  INTEGRATION  OF  DGF/DT   1 

F  SCALAR  FUNCTION  Q  IN  THE  ORIGINAL  DERIVATION  IS  DENOTFD  BY 
IS  PROBLEM  AND  IS  DENOTED  BY  Y(6)  AND  YI7)  RESP.,  IN 
CORAM 

6)  =  (R1*S)/T«-(K1/T)*Y(7)*Y1(2,N)*(1.-YI(2,N)/AN)+Y(7) 
7)=(R2*S)/T+(R2/T)*Y(7)*Y1(2,N)*(1.-Y1(2,N)/AN)-Y(7)*(C+AT(N)  I 
-2.*Y1(2,N)/AN)  I 

N  1 

SFCOND  PART  ENDS 


END 


ee        subroutine  outp(x,y,dery, ihlf,ndim,prmt)  131 

8<5  DIMENSION  PRMT(l0)tCERY(10)fAUX(8fl2)tY0(2tl01),Yl(9,101)fY(l())f.-.r 

1(101  )  ,ATNEW(  101) 


<5C  COMMCN  Yl,AT,ATNEW,YC,AtB,C,AN,F,CIfPI,CA,KfNK,EP,»U,R2,S,l,J, 

SI  IFIK.EC.2)  GO  TO  2C 

C   THIS  PART  OF  THE  SUBROUTINE  IS  FOR  THE  FORWARD  INTEGRATION  ONLY 


S2  IF  (X.NE.O)  GO  TO  21 

<52  J  =  0 

<94       23  FORMAT  (•-   NO.    GRIC.PT.       INVENTORY  SALES 


1      PROFIT  ADVT.') 

95  PRIM  23 

S6       21  J=J+1 


C   STERING  THE  VALUES  OF  STATE  VARIABLES  AT  EACH  GRID  POINT,! 

C   USED  IN  THE  SE  CND  PART  OF  THIS  SUBROUTINE  FOR  THE  CALCULATION  OF 
C   THE  NEW  VALUES  OF  THE  CONTROL  VARIABLE  AT  1C1  GRID  POINTS 


97  DO  22  M=l,2 

96       22    Yi(M,J)=Y(M) 
99  ABC=-Y(3) 


ICC       24  FORMAT  (IH  , I  A, AX , F6. 2 , 5X, A ( E 1 2 . A , 7X )  ) 
1C1  PRINT  2A, J,X,Y( 1), Y ( 2 ) , ABC , AT U ) 

1C2  IF  (J. EC. 101)  PRMT(5)=1. 


1C3  RETURN 

C  FIRST  PART  FOR  FORWARD  INTEGRATION  ENDS 

C 


C   SECOND  PART  FOR  BACKWARD  INTEGRATION  ONLY    

1C4       20  IF  IX.NE.l )  GO  TO  25 

1C5      306  FORMAT  {'-NO.   GRID  PT.     Zl  Z2  Pll  PI 


12  P22  QF1  GF2  S  I' 

2) 
1C6  PRINT  306 


1C7  N=NK+l 

1C6       2  5  N=N-1 

1C9  R1=Y(4)*Y1 (2,N)*( 1 .-Y 1 ( 2 , N ) / AN ) 


IK  R2=-2.*CA*AT(N)+Y(5)*Yl(2,N)*(  l.-Yl(2,N)/AN)+Y(2)*(l.-2.*Yl(?,\)/A 

IN) 
111  S=-2.*CA*Y1 (2,N)*AT(N)+Y(2)*Y1(2,N)*( 1 .-Y 1 ( 2 , N ) /AN ) 


112  T=-2.*CA*Y1(2,N) 

C   BACKWARD  INTEGRATION  OF  DZ/DT   2  EQUATIONS   

C   CALCULATION  OF  THE  NEW  VALUE  OF  THE  CONTROL  VARIABLE  AT 


C  N  TF  GRID  POINT  - 
112  3C  XXI=Y1(1,N)-Y0( I ,N) 
114  XX2=Yl(2,N)-Y0(2,N) 


115  ATNEW (N)=AT(N)-(EP/T)*(S+Y(7)*Yl(2,N)*(l.-Yl(2,N)/AN))-(Rl*XXl+Ri?* 
1XX2J/T 

116  308  FORMAT  ( IH  , I  A, LX , F 5. 2 , IX, 9 ( E12. A f  IX )  )   

117  PRINT  308,  N,  X,  (Y(  II*),  IM=1,7),  S,T 
11£  IF  (N.EO.l)  PRMT(5)=1. 
US      310  RETURN 


C   SECOND  PART  FOR  THE  BACKWARD  INTEGRATION  ENDS 


12C  END 


c 

c 


SUBMri.  [  INT  RKGS(PRMT,Y,DERY,NDIM,  I  HLF, FCT , OUTP , AUX ) 


132 


CIMENSIDN 
X=P«^T( 1 ) 
H=PRNT (  n 


Y  (  1  )f  DERY(  1),  AUX(  R,l),A(4),BU),C(<»),PRMT(l) 


PRMT (5)=P. 

CALL  FCT (X,Y,DERY) 


U 


c 
c 


PREPARATIONS 
A{  l)  =  .5 


FOR  RUNGE-KUTTA  METHOD 


A(2)=. 292893? 
M i)=l. 707107 
A(  4)=.  1666667 


B(  l)  =  2, 
B(2)  =  l, 

!  (  i)  =  l, 


B  (  '♦ )  =  2  . 
C( l)=.5 
C(2I=. 2928932 


C(3)=l. 707107 
C(4)=.5 


l_ 


P-tEPARATIONS  OF 
DO  3  I=l,NDIM 
AUX( 1,1 )=Y(I  ) 


FI'tST  RUNGE-KUTTA  STEP 


AUX(2, I  )  =  DERY(  I ) 
AUX( 3,1 )=0. 
i    AUX(f , I  )  =  0. 


C 

c 


RECORDING  OF  INITIAL  VALUES  OF  THIS  STEP 
7  CALL  CUTP( X,Y, DERY, IR EC , NDI M, PRMT ) 


C 

c 


IF{P"MT( 5) )40,8,40 


START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
8  J=L 
10  AJ=A(J) 


pj=e( j) 

CJ=C( J) 

DO  11  I=l,NDIM 


L 


Ri=E*DEKY( I ) 
R2=AJ*(R1-BJ*AUX(6, I ) ) 
Y(  I  )=Y(  I  )+R2 


R2=R2+R2+R2 
II  AUX(6, I  )=AUX(6, I  )*R2-CJ«"U 
IFU-4)  12,  15,15        


L_ 


12  J  =  J ♦  I 

IF( J-3)  13,  14,  13 


15S 
16C 
It  1 


13  X=X+.5*H 

14  CALL  FCT(Xf Y,DERY) 
GOTO  10 


133 


C 
C 
C 


END  CF  INNERMOST  RUNGE-KUTTA  LOOP 


15    DO    2<3     I=1,NDIM 
AUX{ 1, I )=Y( I ) 
AUX( 2,1 )=DERY( I ) 


29    AUX(6, I )=AUX(3, I ) 

CALL    CUTP(X,Y,DERY,  IHLF  ,  NIC  I  M  ,  PRMT  ) 
IF(PRMT(5)  )40,30,40 


30  DO  31  I=1,NDIM 
Y( I )=AUX(1, T ) 

31  DERY( I )=AUX(2,I ) 


171 
172 
172 


GO  Tf  f- 
40  RETURN 
END 


36 
AC 


Al 


A2 


A3 

AA 


A5 
A6 


A7 
A6 


,  (AT( I ) , I  =  l,NK) 
DO  515  JK«1 , 101 

515  y    ,   )  -Y  (  h,  JKJ  /ICO. 

Arc  FOLI       VALUES  OF 

READ    IN' ) 
1      6  J 
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6  *    l 

1     y/AR 

i   !  <r* 

P  <  I  v 
LT2(  1 


THE    6    STAff"     WMAuLFS     AND    3    CONTRf'L 


61 


1     63  

/H     (IH    ,  (  IH    , I3,1X,9(F10.5,IX) )/) 

I    61,(1 ,Y0( LtI )tY0(2, I ),Y0(3,I ),YO(A,I ),Y0<5, I) ,Y0(6, I), 

)  ,AT( I ) ,I  =  1,NK) 


T  1  (  I 


nn  loop  for  i r-  at  ions 


MAIN 

on  ico  u=ifrrMAx 

AI   (•   ************************************ 


iTFP.a  if  •• 


1     '.lA,'     ***********************************») 

[NT  200 i U 

PRMT ( 1)=0.0 


AS 
51 


PKMT 
PRMT 

I 


?)  =  l. 

3  )  =  .  0  1 
A)  =  .l 


52 
*>3 
5A 


NDI^=7 

DERY ( I )=MD IM 

DC^Y( I  )=1.E0/DERY(  1 ) 


55 
56 

S7 


DO  I 

I  DEKY 

Yd) 


l=?,NDIM 
(  I  )  =  DFKY(  1) 


5S 
6C 


Y(2) 
Y(  3) 
Y(A) 


.A3 
.5  I 

.A3 


61 

Y(^)=R. 

62 

Y  (  6  )  =  .  0  1 

6  3 

Y( 7)=0 

6A 

KSL=1 

65 

Al 

CALL  «KGS(PRMT,Y,DERY,NDIM,  IHLF  ,  FC  I"  ,OUTP  ,  AUX  ) 

66 

PRM  (  1  )  =  l  . 

67 

T(2)=0. 

|~ 

66 

PRMT(  n^-.oi 

6S 

PRMT(A)=10. 

7C 

NDIM=33 

71 

:Y(  1  JsNDIM 

72 

0ERY(1)=1.E0/DERY(1) 

73 

DO  3  l=?,NDIM 

7A 

3 

DERY  (  I  )=UFRY  (  1  ) 

75 

DO  A  1=1,33 

76 

A 

Y(  I)  =  0. 

77 

KSL=? 

76 

A6 

CALL  Rt  ,,( PRMT, Y, DERY, NDIM, IHLF,FCT,OUTP,AUX) 

7S 

DG  120  L=l ,NK 

ec 

Y'U  1,L)=Y1  (  1,L) 

81 

l2tL )=Y1 (2,L) 

82 

YU( 3,L )=Y1  (  1,L) 

83 

YMA,L)=Y1(A,L) 



L 

84 

85 
86 


Y0(5fU=Yl (5, L) 

YIM6,L)=YI  (6,L) 
r  KL)  =  TINEW(L) 


87  T2(L)=T2NEW(L) 

£8      120  Ar(L)=ATNEW(L) 

8S       100  C0M1 IMUh  


9C 
91 


STOP 
END 
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92 

93 


SUBROUTINE    FCT    (X,Y,DERY) 

,  RMT(10),Y(42)  ,DERY(42),AUX(8,43),YU40,202),tWeW(202) 

II  (  I ,ATNEW(2 u2) ,T1 (202) ,  T2(202) , AT ( 202 ) , YU ( 6, 202 ) , A( 10) tL( 1-^ 

2),M(10)  J 


2),N(1C) 

S4  COMMON    YlfTl  »T2NEWtATNEWvTlvT2f AT9NKtXXltXX2f XX3vXX4fXX5»XX6» Y 

L,FTl,FT2tFT3tFT4,FT5>Sli        i      jj   '  L  ,  J/  t  ■'  ),R4 ,  R5  ,  R6 1  R7,  R8  ,  R9  t  RlU  ♦  R  1 1  *- 
;i/,KM,Kl4,R15,Kl6,Rl7,Rl8fA,J,NfEP,KSL,R,DIS,Vl,\/2,CQ,CfAN,AIN<fT 
»CltC2tC3fCIffCA«CTtl    •.       iXI , YI « GA,GB,RTR1, RTR2f RTR3f RTR4,R1      >tRT 
<>6tRTR7>RTR8tRTR9,RTRl      ,  1 L ,RTRI2>RTRI3>RTRI4tRTR15,RTR16tRr«17    J 

9  5  MCN    RTR18VRTR19,RTR20VR1  12  Li  U  r>  1  ,  i:  TS2  ,  RTS3,  RTS4,  RTS5,  R  rS6,:<TFU~ 

LfRTFC2tR1  tTFQ4fRTFQ5,RTFQ6,FQl,FQ2tFQ3,FQ4tFQ5,FQ6,  ■?  , 

1     II,'    '  r2,KMl,BM2,BM3tBM4tBM5t  GlrG2tG3tG4 - 

[I      IKSL-EQ.2I    GO    TO    2  5 
97  [f     (X.NE.OJ    GO    TO    42 

Sfc J=l - 

9S  A 2    DERY(1)=(DIS/V1 )*(XI-Y(i) )-GA*EXP (-EA/ ( R * T L ( J)  )  )*Y( 1) 

ICC  DERY(2)=<DIS/V1)*<YI-Y{2))-GB*EXP(-EB/«R*T1(J)))*Y(2)+GA*EXP(-EA/I 

LR*T1(J) ) )*Y( 1) 


DERY(3)=(DIS/V2)*(Y(1)-Y(J))-GA*EXP(-EA/(R*T2(J)))*Y(3) 

UE*Y(4)  =  (DIS/V2)*(Y(2)-YK)  )  -GB*EXP  (  -EB/  (R*T2  ( J )  )  )*Y  (4  )+GA*EXP(- 

1  /(R*T2(  J)  )  )*Y(3) - 

(S)=hIS*Y(4)-CQ*Y(6) 
DERY(6)=(C+AT( J) )*(Y(6)-( (Y(6)**2)/AN) ) 

DERY( 7)=CQ*C1*Y(<S)+L2*DIS*Y( 3 ) +C3*U I S* ( 1-Y (  3 ) -Y ( 4 ) ) -C I  * (  (AIM-YC5)  }- 
1**2 ) -C A* ( AH J)**2)*( Y(6)**2)-CT*( ( T LM- T 1 ( J ) ) **2+ ( TL ( J ) -T 2 ( J ) ) **2 ) 
RETURN 
2b     IE  (  X.NE.l  )  GO  TO  W 


N=NK 
*7  CALL  CALCL (X,Y, DERY  ) 
INTEGRATION  OF  DZ/DT 


CF«Y( 1  )=-BMl*Y( 1  )-Y(2)*GA*EATl-(DIS/V2)*Y(  3) 
Y(2)=-RW2*Y(2)-Y(4)*(0IS/V2) 
'Y(  U=DIS*(C3-C2)-HM3*Y(  3  )-Y  (  4  )  *  G  A*E  AT  2 


DERY<4)=C3«DIS-Y(4)*BM4-DIS*Y(5) 
DERY(5)=2.*CI*(YK5,N)-AIM) 

DERY  (  6  )  =-C  1  *C0+2  .  *C  A* ( A  T  (  N  )  **2  )  *  Y  1  (  6,  N  )  +Y  (  6  )  *CQ*(C  +  AT  (N)  )  *l  l.-2.*V 


LK6,N)/AN) 

INTEGRATION  OF   DP/DT     

UERY(7)=-2.*Y(7)*BM1-2.*Y(8)*GA*EAT1-2.*Y(9)*(DI S/V2)+RTH1 


DERY(R)=-Y(8)*(BM1+BNV)-Y(  I  3  )  *GA*b  AT  I  -  (  Y  (  I 't  ) +Y  (  1  0  )  )  *  (  D  I  S/ V/' ) +RTR  2 
DERY(9)=-Y(9)*(BM1+BM3)-Y< 10) *GA*EA T2-Y ( L4)*GA*EAT1-YU8 )*(0I S/V2) 
UR1 


DERYI  10)=-Y(  10)*(BMI+BM4)-Y(15)*GA*EAT1-Y(  19)*(UIS/V2)-Y(  11  )*DIS+R 
1TR4 
DERY(  ll)=-Y(ll)*BMl-Y( !6)*GA*EAri-Y(20)*(DIS/V2)+RrRS  


DERYI  12)=Y( 1 1 )*CQ-Y( 1 2 ) * ( BM 1 +BM5 ) - Y ( 1 7) *GA*C AT  1- Y ( 2 1 ) * ( D I S/V2 )  ♦  R T P 
16 

(YllO=-2.*Y(  13)*BM2-2.*Y( 1 5 ) * ( D I S/ V2 ) +RTR7 

■Y(l^)=-Y(1A)*(BM2*BM3)-Y( 15 ) *GA*EAT2-Y ( 19)*(DIS/V2)+RTR8 
DERY(  15)=-Y( 15)*(BM2  +  BNA)-Y (22 )*(UIS/V2)-Y(  16)*0IS  +  RTR  * 


125 
126 
127 


12E 
129 
13C 


:Y(  16  )=-Y(  16)  *BM2-Y(  23  )*(DIS/V2)+RT, UO 
DERY(17)=Y ( 16)*CQ-Y (17)*(BM2+BM5)-Y(24)*(DIS/72)+Rl,l! 
DERY(18)=-2.*Y< 1 8 ) *BM 3-2 . *Y ( 19)*GA*EAT2+RTR12 


138 


DERY(  19)=-Y{ 19)*(BM3+BM4)-Y(22)*GA*EAT2-Y(20)*DlS-MTRl3 

DERY(20)=-Y(20)*By3-Y(23)*GA*EAT2+RTR14 

nEKY(21)=Y(20)*CQ-Y(21)*(BMi+BM5)-Y(?4)*GA*EAT2*RTRlS 


131 
132 
133 


DERY(22)=-2.*Y(22)*BM4-2.*Y( 23)*DIS+RTR16 
DERY(23»=-Y(23>*BM4-Y(25)*DIS+RTR  I  1 

DERY(24)=Y(23)*CQ-Y(24)*(BM4+BM5)-Y(26)*DIS+RTR18 


134 
135 
136 


DERY(25)=RTR19+2.*CI 
DERY(26)=Y(25)*CQ-Y(26)*BM5+RTR20 

DERY(27)=2.*Y(6)*(C+AT(N) )/ AN+2. *CA*( AT < N ) **2 ) +2 . *Y ( 26 ) *CQ-2. *Y< 2 


137 


l)*BN5+RTR2l 


INTEGRATION    OF    DQ/DT 

DERY(28)=RTS1+RTFQ1-FQ1 


136 
13S 
14C 


0ERY(29)=RTS2+RTFQ2-FQ2 
D£RY(30)=RTS3+RTFQ3-FQ3 
DERY(3l )=RTS4+RTFQ4-FQ4 


141 
142 
143 
144 


DERY« 32 

OERY{ 3  3 

RETURN 
END 


)=RTS5+RTFQ5-FQ5 
)=RTS6+RTFQ6-FQ6 


su 

(I 

1  I  2 


UTINI 

: 

(202) 


14  / 

L.FT1, 

^,Ll  , 

A6,  '  1  R 

146 

4MC 

1  .  MFC 

i  r  i  v  e 

14S 

ir  (k 

15C 

II  (X 

151 

J  =  0 

152 

U  J  =  J* l 

DOTiM  X,  Y 

PRMT  (1      )  , 

,ATNEH 120 

,  T 2N 

r  3  «  F  T  4  9  F  T 

,  '1>.  U6,l 

ItCIfCAtCT 

,  tTR9,  tTI 

18,RTR1  )_i 

-.  >,  <TFQ4, 

Ml , HM2.BM 


V  ( 42  1 

2  )  ,  T  1 
[W.AF 
5,i>l, 
R 1 7  ,  R 


,  IHLFi   Ji)IM,PKMT) 

'Y('.2),AUX(8, 

(2  02 ) ,  \2  (202) ,AT 


«  n  iYiC40*202) ,  ri?8  .  i 

2 ) ,Y0(6, 2  02), A(  lu) 


»T1, 

t  S3  *  R 

18, A, J, 


.2) 

I    60   Tl       I 


tEA,E 
K  1 U  ,  R 
RTR20 

3,  BM4 
36 


BfXI«Yl 

,RTR21, 


T2,AT 
LtR2, 

NtEP, 


,RTI 

,BM5,G1 


,GA,G 
R12,R 
RTSl, 


) , AT (202) ,Y0(6,202), A(  lu) 
,NK,XX1,XX2,XX3,XX4,XX5,XX6, 
, (4vR5ffR6tR7, R8, R9, R10.R11 
KSL,",MS,V1,V2,CQ,C,AN, A  I  . 
,  ITR1 »RTR2tRTR3vRTR4vRTR5,R 
TR13,RTR14,RTR15,RTR16,RTR17 
RTS2,RTS3, RTS4.RTS5.KTS6.RTF 


) 
1 


FQl.F 
,G2,G 


Q2,F 
3,G4 


R1,RTR2,RT  13VRTR4«  tl  15  ,R 
,RTR14,RTRi5,RTRl6,RT  tl7 

tRTS3,RTS4,RTS5,RTS6,.<TF', 
Q3,FQ4,FQ5,FQ6,FAT1,EA1 


1 


32 
15 


DO    52 

Yl(Kr 


K=l,6 
J)=Y(K  ) 
f     (1H    ,  I3,3X,F6.3,2X,7(C12.4,  iX) ) 


16  PRINT     15, J,X,  ( Y(  I  )  ,  1  =  1, 7) 
IF     U.EQ.NK1     PRMT(5)  =  1. 

17  RFTURN 


36 


IF     (  x 

!K  + 
33    N=  ,-1 


.NE.l )     GO    TO 

I 


33 


CALCL(X,Y,DERY 

T     (  1  H     ,  15  ,  L  X , F 

20,N,X, (Y( I ) , 


L 


20 


CALL 

FORM 

PRINT 


20fN,X, (Y( I), 

20, N, X, (Y( I ) , 
20,N,X,  (Yd), 


) 

5. 5, 1X,8(E12 

1=1,8) 


.4,  3X)  ) 


P  UNT 
PRINT 
PRINT 


1=9, 16) 

1=17,^4) 

1=25,32) 


(DF'/DTHETA). (Q) 


7T 


91 


DO    91 
Yl(  14 


1=1,33 

7,N)=Y( I ) 


75    FOUl  =  FTi*Yl(  35,N)  + 
FD1 C2=FT3*Y1 (37,N)+ 
FOTQ «  =  FT5*Y1U0,N) 


FT2*Yi(36fN) 
FT4*Y1 (38, N) 


S+(DF'/DTHrTA).Q 


SFDla 

SF02  = 


Sl+FDTQl 
S2+FDT02 


18L 
181 
182 
183 


SF03=S3+FDT03 


T(  INVERSE 

TISFCL=A( 1 )*SFD1+A(4)*SFD2+A(7) 


)  .  (b+(OF«/OTHETA) .Q) 
*SFI,3 


I  ISFC 
TlbFC 
f-P  =  .l 


2=A(2)*SF01+A( 
3=A( 3)*SFD1+A( 


5)*SF02+A( 8)*SF03 
6)*SFD2+A(9)*SFb3 


TIKU 

T  I  <4  = 

r  i 

I  i  <  i  c 


T( INVERSE 

A( i)*Rl+A(4)*K2*A( 7)*R3 
A(  i  )*R4+A(4)*K5+A(7)*R6 


).R 


A(  l)*R7  +  A( A)*R8  +  A( /)*K( 

=  A(1)*R10  +  A(4)*R1H-A(7)*RI, 


184 
185 
186 


TI-U3  =  A(l)**13  +  AU)*R14«-M7)*R15 
riU6=A(l)*R16+A(4)*Rl7+A(7)*R18 
TIR2=A(2)*R1+A(5)*R2+A(8)*R3 


1U0 


187 

i  e  e 

189 


19C 

IS  1 

192 


T IR5=A(2)*R4+A(5)*R5+AC8)*R6 

TI  18=AC2)*R7+A(5)*R6+A(8)* 

1  Ull  =  A(2)*R10+AC  5)*Rll+AC8)*Rl2 


TIR14=A(2)*R13+A( 5)*R14+A(8)*R15 
TI  U7  =  A(?)*R16  +  AC  5)*R17+AC8)*R13 
TIR3=AC3)*R1+A{6)*R2+AC9)*R3 


193 
194 
195 


riR6=Am*R4  +  At6)*«5  +  AC9)*R6 
TIR9=A(})*R7+AC6)*R8+A(9)*R9 
r  Ml  2  =  AC  3)*R10  +  A(6)*RiHvA(9)*Rl2 


196 
197 


TIR15=AC3)*R13+A(6)*R14+A(9)*R15 
TIR1P=A(3)*R16+A(6)*R17+A(9)*R18 
X  I J  +  l  )-X ( J  ) 


196 
199 
2CC 


XXl=Yl ( l,N)-YO( i, N) 
XX2=Yi(2,N)-YOC2,N) 
XX3=Y1( 3,N)-YQ(3,N) 


2C1 
2C2 
2C3 


XX4=Y1(4,N)-Y014,N) 
XX3  =  Y1(5,NJ)-Y0C5,N) 
X  X6=Y 1 ( 6,N)-Y0(6iN) 


2C4 

2C5 
2C6 


r  UX1  =  XX1*TIR1  +  XX2* 
TIRX2=XXI*TIR2+XX2* 
TI  vX}=XXl*TIR3+XX2* 


TIR4+XX3*TIR7+XX4*TIR10+XX5*TIR13+XX6*IIR16 
TIRS  +  XX3*TIR8  +  XX4*TI!U1  +  XX5*TIR14  +  XX6*TIR17 
TIR6+XX3*TIR9+XX4*TIR12+XX5*TIR15+XX6*TI 


2C7 
2C8 


TINEWCN)=T1 (N)-EP*T 
rZHEM CN)=T2 (N)-EP*T 


IMPROVED  VALUES 
ISFQ1-TIRX1 
1SF02-TIRX2 


OF  CONTROL  VARIABLES 


2C9 
21C 
211 


ATNEW  (N)=AT(N)-EP*f 

21  FORMAT  (1H  ,I3,3X,F 

PRINT  21,N,Xf YC33), 


ISF03-TIRX3 
5.2,2X,4CE12.4,2X) > 
T1NEMN)  tT2NEW(N)tATNEW(N) 


212 

213 
214 


19 


IF  (N 

RETLR 


.EO.l  ) 

N 


PRMTC5)=L. 


2  15 
216 


217 


218 


219 
22C 
2/1 


222 
223 
224 


22* 
22t 
221 


DI 
11  , 
2)  , 


rui  I 'si 
NSION 
v  1202 


\Li.l.  (X,Y,  DERY  )  lUl 

4T(l0),Y(42),DERY(42),AUX<8,43)f Y1(40,202),T1I 

>,UNFW(2u2),ri(2(v2>,T2<2c2)tAT<2C2),Y0(6,2C2),/WlC) 


,L(  It 


10) 
CM  Yl 
I ,i I  t»FT2, 
1  3,Rl 
,C2.C 
R7,RT 


2 12 1 
3M,C1 


,  T  1  NEW, 

r  i  >,FT4 
4,  tl5vR 

3, CI  »CA 
,RTR9 


T2Nt 
,FT5 


X6,Yf 
R 1 1  , 1 


Chmp 
l.KIF 

111, 


ton 


rN  RT 
C2.RT 

=  PXP( 
=  FXP( 
=  EXP( 


R1R,RTR 

FQ3.RTF 

BM1, 

-  F  A  /  (  I  * 

-FA/ (R* 

-EB/(R* 


I  . 
BM2  = 


=  FXP( 
-(DIS 
-(OIS 


-EB/(R* 

/  v  i )  - G  a 

/Vl )-GB 


T2(N 
*EAT 
*EBT 


)  )  ) 
1 

1 


BM5= 


-(CIS 

-  (IMS 
I  C  +  M 


/V2)-3A 
/V2 )-GB 
(N)  )*( 1 


►EAT 

*EBT 
.-2. 


2 

2 

*Yi (6,M)/AN) 


226 
22S 
23C 


51=  < 

G2=( 


-CA/( 
-FB/( 
-E/\/( 


R*(  TKM 
R  *  C  T 1  ( N 

R  *  (  T  2  ( "i 


)  **2 
)**2 
)**2 


)  )  )* 
)  )  )* 

)  )  )* 


GA*3ATl 
GB*EBT1 
GA*FAT2 


231 
232 
211 


G4=(-EB/(R*(T2(N)**2)) )*GB*CBT2 
FTI=Y1 ( 1,N)*G1 
Fr2=Yl(2,N)*G2-FTl 


234 
235 
236 


FT  3= 
F  T4>= 
FT5- 


Yl(  3, 
Yl(4, 
Yl(6i 


N)*G3 

N)*G4-F 
N)*(l.- 


T3 
Y1(6,N)/AN) 


237 


ALL  MATRICES  TREATED  AS  VECTORS  ,  CfLUMNWlSE 

MATRIX   R   (3X6)   

R1=Y(7)*FTI+Y(8)*FT2+Y(1)*G1-Y(2)*G1 


236 
2  3S 
24C 


R2  =  Y 
R  3=Y 
R4  =  Y 


(9)*F 
(12)* 
(8)*F 


T3  +  Y(  10 
FT5 

ri+Y( 13 


)*FT4 
)*FT2+Y(2)*G2 


241 
242 
243 


R5  =  Y 
R6=Y 
R7  =  Y 


(  14)* 
(17)* 

m*F 


FT3+Y( 1 

FT5 

Tl  +  Y(  14 


5)*FT4 
)*FT2 


244 
245 
246 


R8=Y 

R  >  =  Y 

'  !   = 


(IB)* 
(21)* 
V(iOJ 


F  T3  +  Y( 1 

FT5 

*FTl+Y( 


9)*FT4-G3*Y(4)+G3*Y( 3) 
15)*FT2 


247 
246 
24S 


25C 
251 
252 


253 
254 


Rlis 

R12^ 
Rl  3^ 


Y(IS) 
Y(24) 
Y(  11) 


*FT3+Y( 

*FT5 

*FTl+Y( 


22)*FT4+G4*Y(4) 
16)*FT2 


R14: 
RIO: 


Y(20)*FT3+Y(23)*FT4 

Y(26)*FT5 

Y(  12)*FT1+Y( 1 H*FT2 


*ir 


Y(21  )*FT3  +  Y( 24)*FT4 

Y(27)*FT5  +  Y(6)*(l.-2.*Yl(6,N)/AN)-4.*CA*AT(N)*Yl(<s,, 


25 


Si=2.*CT*( 
12*Y(2)-Y1 ( 


MATRIX     S   (  1X1)  

riM-TKN)  )+2.*CT«(TKN)-T2(N))*Yl(ltN)*Gl*YI 
L,N)*G1*Y(2) 


£1*2 

1) +f L( 2 


. 


256 
257 


258 


261 
262 


265 
266 


S2=2.*CT*( 
1*Y{4) 
S3=-2.*CA* 


T1(N)-T2(N))*Y1(3,N)*G3*YM)+YL ( 4 , N ) *G4*Y ( 4 ) - Y 1 (  , 
AT (N)*( Yl  (6,N)**2 )+Y(6) *( Yl ( 6t N )- ( Y 1 ( 6 , N ) **2 ) / AN ) 


DT1=-2.*(Y 
U-2. »Y1(2, 


MATRIX     T   (3X4) 

1  (  l,N)/Tl  (N)  )*G1*Y{ l)  +  ( Y(  1  )*G1* 
N)*G2*Y(2)/Tl (N)+EB*Y1 (2,M)*G2* 


EA*Y1 (l,N))/( 
Y  (  2  )  /  ( R  *  (  F  1  (  i 


R*(Tl(N 

)**2  I  )  + 


)*<V  ) 
(2.*' 


21  (  l,N)*Gl*Y(2)  )/Tl(N)-FA*Yl U tN)*Gl*Y (2 ) / (R*( Tl (N)**2) ) 

259  DT2=-2.*CT 

26C  DT3=C.  


DT4=2.*CT 
UT5=-2.*CT 
L2)  )-?.*Yl( 


-2.*Y1(3,N)*G3*Y(3)/T2(N)+FA*YI 
4»N)*G4*YU)/T2tN)+EB*Yl  U,N)*G 


( 3,N)*G3*Y( 3) 
4*Y(4)/  (R*(T2 


/«R*(T2 

(  1)**2  ) 


I  I 
»+2. 


2YH3iN)*G3«Y(4)/T2(N)-EA*Yl(3iN>*Y(4>*G3/(R*<T2«N)**2)) 

262  DT6=0. 

264         or7=c.  


DT8=C. 

0T9=-2.*CA*( Yl(6tN)**2) 


MATRIX   R'T 


(6X3) 


267 
266 
269 


RT1=R1*DTI+R2*DT2+R3*DT3 
RT2  =  R4*DTl  +  R5*DT2  +  r,v6*ur3 
Rr3=R7*0Tl+R8*DT2+R9*DT3 


27C 
271 
272 


RT4=Rlu*nTl+Rli*DT2+R12*DT3 
RT5=R13*DT1+R14*DT2+R15*DT3 
Rro=Rl6*DTl+R17*0T2+R18*0T3 


273 
274 
275 


RT7  =  R1*DT4  +  R2*DJ"5  +  R3*DT6 
RT8=R4*0T4+R5*QT5+R6*DT6 
RT9=R7*DT4+R8*0T5+K9*DT6 


276 
277 
276 


RT10=R10*DT4+Rll*DT5+R12*DT6 
RT11=R13*DT4+R14*DT5+R15*DT6 
Rri2=R16*0TA+R17*Dr5+R18*DT6 


279 
28C 
281 


RT13  = 
RTI^h 
RT15' 


R1*DT 
R4*DT 
R7*DT 


7  +  R2 
7  +  R5 
7  +  R8 


*DT8 
*DT8 
*DT8 


+R3*0I9 
*R6*DT9 

+R9*DI  ) 


282 
283 
284 


RT16=R10*DT7+R11*DT8+R12*DT9 
RT17=R13*DT7+R14*DT8+R15*DT9 
RT18=R16*DT7+R17*DT8+R18*DT9 


285 


286 
287 

268 


289 
29C 
291 


292 
293 


M  MR  I 
RTR1 


X       DP/ 
=R1*RT 


RTR2= 

RTR3; 


Ri*RT 
R1*RT 

R1*RT 


DT  E 
1+R2 
2+R2 

i  +  R2 
4  +  R2 


MATRIX   R'TR 
EING  SYMMETRICAL  ONLY 
*RT7+R3*RT13 


(6X6)      

HALF  THE  MATRIX  R'TR  lb  fAKEN 


*RT8 
*RT9 

*^T1 


+R3*Rr 14 
♦R J*RT15 

0+R3*RT16 


I  I     '- 
RTR6^ 
RfR7^ 


R1*RT 
R1*RT 
R4*RT 


5  +  R2 

6  +  R2 
2+R5 


*RT1 
*RT1 

*RT8 


1+R3*RT17 
2+R3*RTl8 
+R6*RT14 


RTR8^ 
RTR9^ 


R^*RT 
R4*RT 


3  +  R5 

4+R5 


*RT9 
*RT1 


+R6*RT15 
0+R6*RTl6 


254  'i     i         t4*RT5-*-R5*RTlH-R6*RT17  II13 

255  RTKl  l=R4*RT6*R5*Rr  1?*K6*KT18 

I   'l?=R7*M  1  >+R8*R  f94R  >*R  I  I '» 


2S7  'I  f*kT'»4Ra*RTl(J*R9*RTl6 

2<5E  f*RT5+R8*RTll*R9*RTl7 

2SS  LI    :  lr=R7*RT6+R3*RTl2+R9*RT  [M 


3CC  10*RT4  +  'U  L* R T  10*R  12*RT 16 

501  I    '  1  7  =  R10*RTC>  +  RL  l*!<  ri  L+R12*RT1  7 

3C2  RrRiP=RiO*RT6*Rll*RT12*Kl2*RTi8 


3C3  !U9  =  Rl3*RT,i+Rl4*RTlH-Ri5*RrW 

3C4  Rr,R2C  =  R13*«T6  +  ^l^*RTl2  +  '<l,>*RT18 

i-Cc  i     /  1  -  ;16*RT6*R17*RT12*.n1«*RT18 


C CALCULATION    OF        T-INVERbF: 

3ce  Ml)=DTl 

3C7  A(2)=DT2 


j oe  a  (  n  =  d  t  3 

3CC  A(4)=DT4 

3  1  C  A  (  ',  )  =  Q  T  S 

311  MG)=DT6 

3  12  A(  7) =01 7 

313  A{8)=DT8 


-iWr- 


341 

RTF1=RTI  l*FTl 

342 

RTF2=RTI2*FT1 

343 

RTF3=RT13*FTl 

344 

RrF4=RTI4*FTl 

345 

RTF5=RTI5*FT1 

346 

RTF6=RTI6*FT1 

347 

RrF7=Rn  i*f  r? 

346 

RTF8=RT12*FT2 

34S 

UF9  =  RTI3*FT2 

35C 

RrFlC=Rri4*FT2 

351 

RTFli=RTl5*FT2 

352 

RrF12=RTI6*FF2 

353 

ri.rF13=«ri7*FT3 

354 

RTF14=RTt8*FT3 

35  5 

RTF15=RTI9*FT3 

356 

^TF16=RTI10*FT3 

357 

RTF17=RTU1*FT3 

356 

RTF18=RTI12*FT3 

35<; 

RTF19=RTI7*FT4 

36C 

RTF2C=RTI8*FT4 

361 

KrF2l=RTl9*FT4 

362 

RTF22=RTT10*FT4 

363 

RFF23=RTlll*FT4 

364 

RTF24=RTI12*FT4 

365 

RTF25=0. 

366 

RFF26=0. 

367 

RFF27=0. 

366 

RTF2P=0. 

36S 

RTF29=0. 

37C  RTF3C=0. 

371  RrF31=RTI 13*FT5 

372  RrF32=RTI14*FT5 


373  RrF33=RTll5*FT5 

374  RTF34  =  RU  16*FT5 

375  RTF35=RTU7*FT5 


376  RTF36  =  RTU8*FT5 

c       

377  RrFCl  =  Y(28)*RTFl  +  Y(  29 ) *RTF 7+Y { 30 ) *RTF 1 3+ Y (  J  1  )  *,<TF  1  9  +  Y  (  3?  )*RTF25+Y 


L33)*RTF31 

376  RrFQ?=Y(28)*RTF2+Y(29)*RTF8+Y( 30)*RTF14+Y< 31)*RTF20+Y1  J2)*RTF26+Y 

I  1  3)*RTF32 


37S  RTFQ3=Y(28)*RTF3+Y(29)*RTF9+Y(30)*RTF15+Y(31)*RTF2H-Y(32)*RrF27*Y 

133)*RTF33 
38C  RTFC4=Y(28)*RTF4+Y(29)*RTFl0+Y(30)*RrF16+Y(Jl)*RrF2  2+Y(32 )*R1 \   28+' 


1 ( 33)*RTF34 
381  RrFC5=Y(2P)*RrF5+Y(2  9)*RTFll+Y(30)*RTF17+Y(31)*RTF23+Y(.V)^TF29+ 
1<  3  3) *  RTF  35 


382  RTFQ6=Y(28)*RTF6+Y(29)*RTF12+Y(30)*RTF13+Y(3l)*RTF2  4+Y(32)*RTT 

1(33)*RTF36 


L    

3e?  )1»Y(28)*BM1 

384  -Y  (.M)*GA*EATl«-Yt2'J)*dM2 


385  i  .  J  Y(,>  OMDIS/V2  )  +Y(  iO)*H(- 13 

386  v (,>9)*<DIS/V2) +Y(  iO)*GA*FAT2*Y(  51 )* 

387  FQ5=Y(3l)*DIS 


388  FQ<   -Y(iJ)MC  +  AT(N))*(l.-2.*Yl(6,N)/AN) 

38<5  RFIuMN 

39C  JD 


C 
C 

c 


C  SUBROUTINE  MINV 

C 

C  IRPOSE 


C  INVERT  A  MATRIX 

C 

C         USAGE 


C  CALL  MINVI A,N,D,L,M) 

C 

C         DESCRIPTION  OF  PARAMET 


C  A  -  INPUT  MATRIX,  DESTROYED  IN  COMPUTATION  AND  REPLACED  \ 

C  R^SULTANI  INVERSE. 

C  N  -  ORDER  OF  MATRIX  A 


C  D  -  RESULTANT  DETERMINANT 

C  L  -  WORK  VECTOR  OF  LENGTH  N 

C  M  -  WORK  VECTOR  OF  LENGTH  N 


C 

C  REMARKS 

C  IRIX  A  MUST  BE  A  GENERAL  MATRIX 


C 

C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  RE^UlRtD 


C  'ME 


C 

C         METHOD 

C THE  STANDARD  GAUSS-JORDAN  METHOD  IS  USED.  THE  DETERMINANT 

C  IS  ALSO  CALCULATED.  A  DETERMINANT  OF  ZERO  INDICATES  THAT 

C  THE  MATRIX  IS  SINGULAR. 

_C 

c  

c 


391  SUBROUTINE  MINV ( A, N,D,L , M )  1U6 

3  92  DIMENSION  A ( 1  ) , L (  1  )  r M  (  1  ) 

c 

c  

c 

C IF  A  DOUBLE  PRECISION  VERSION  OF  THIS  RUUTINE  IS  DESIRED,  THE 

C         C  IN  COLUMN  1  SHOULD  BE  REMOVED  FROM  THE  DOUBLE  PRECISION 

C         STATEMENT  WHICH  FOLLOWS. 

C 


C      DOUBLE  PRECISION  A , D, B I GA, HOLD 
C 

C THE  C  MUST  ALSO  BE  REMOVED  FROM  DOUBLE  PRECISION  STATEMENTS 

C         APPEARING  IN  OTHER  ROUTINES  USED  IN  CONJUNCTION  WITH  THIS 
C         RCUTINE. 

C 

C         THE  DOUBLE  PRECISION  VERSION  OF  THIS  SUBROUTINE  MUSI  ALSO 
C         CCNIAIN  DOUBLE  PRECISION  FORTRAN  FUNCTIONS.   ABS  IN  SI  I 

C IC  MUST  PE  CHANGED  TO  DABS. 

C 

c         

c 


C  SEARCH  FOR  LARGEST  ELEMENT 

C 
393  D=l. C 


394  NK=-\ 

395  DO  8C  K=1,N 

396  NK=NK+N 


397  L(K)=K 

396  M(K)^K 

399  KK=NK+K 


4CC  dI~;A=A(KK) 

4C1  DO  20  J  =  K,N 

402  IZ=N*(J-l) 


403  DO  20  I=K,N 

404  IJ=I7+1 

4C5 10  IF(  ABS(BIGA)-  ABS(A(IJ)))  15,20,20 


406  15  BIGA  =  A(  IJ) 

407  L(K)=I 
406  M(K)=J 


4C9       20  CONTINUE 
C 
C  INTFRCHANGE  ROWS 


41C  J=L(K) 

411 IF(J-K)     35,35,25 


412  25    KI=K-N 

413  DO    30    1=1 fN 

414  KI=KI+N 


415  HOLD=-A(KI) 

416  JI=KI-K+J 


41  1  [)=A(  JI  )  ll|7 

A  (  J  I  ) 

c 

C  INTERCHANGE  CULUf- 

C 

Ms  )  

42C  IF( I-KI  4S,45, 38 

42  1  38  JP=N*(  1-1  ) 

A2  2   DO  4C  J=1,N  


42  2  Jn-NK  +  J 

42  *  JI=JP*J 

425  HOLC=-A(JK) 


426  A  {  J i\  )  =  A  (  J  I  ) 

427  4(,  A(  JI  )  =  HOLU 


C         CIVIDF  COLUMN  BY  MINUS  PIVOT  (VALUE  OF  P  I  VC  I  ELEMEN1  IS 
C         CCNTAINED  IN  *IGA) 


C 


4  26       4  5  IF(BIGA)  48,46,4  8 
425        46  C  =  O.C 
4JC  RETURN 


431  48  DO  55  1=1, N 

432  IF(I-K)  50,55,50 

4  \  i 50  I  ■     1 1 

4*4  A( IK)=A( IK)/(-HI3A) 

435       55  CONTINUE 


C         REDUCE  MATRIX 
C 
4  36 DU  65  I  =  1,N 


437  Its  =  NK  +  I 

436  HOLD=A(IK) 

435  IJ=I-N 


44C  DO  65  J=1,N 

441  IJ=IJ+N 

442 IMI-K)  60,65,60 


443  60  IF(J-K)  62,65,62 

444  62  KJ=IJ-I+K 

445  A(  IJ)=MOLD*A(KJ) +A(  IJ) 


446       65  COHINUE 
C 
C  IVIDE  ROW  BY  PIVOT 


447  KJ=K-N 

44  6 DU  75  J=I,N 


4  45  KJ=KJ+N 

45C  IF(J-K)  70,75,70 

451        70  A(KJ)=A(KJ )/BIGA 


452       75  CONTINUE 


C  PRODUCT  OF  PIVOTS  lU8 

C 
452  D=D*fHGA 


C 

C 

c 

PLACE  PIVOT  BY  RECIPROCAL 

454 

455 

c 

80 

A(KK)=1.0/BIGA 
CONTINUE 

456 

c 

c 

1  TNAL  ROW  AND  COLUMN  INTERCHANGE 

457 

458 
45S 

100 
105 

IE(K)  lr>0,150,L05 
I=L(K) 

46C 
461 
462 

108 

IF(I-K)  120,120,108 
JQ=N*(K-1) 
JR=N*(I-i ) 

463 
464 
465 

DO  110  J=1,N 

JK=JG+J 

HOLD=A(JK) 

466 
467 
466 

110 

Jl=JR+J 

A( JK)=-A( JI  ) 
A(JI)  =HOLD 

46<5 
47C 
471 

120 
125 

J  =  M(K  ) 

IE(J-K)  LOO, 100, 125 

KI=K-N 

472 
473 
474 

DO  130  I=1,N 
KI=KI+N 
HOLD=A(KI ) 

475 

476 
477 

130 

JI=KI-K+J 
AIM  )=-A(  J  I  ) 
A(JI)  =HOLD 

476 
47S 
48C 

150 

GO  TC  100 

<ETURN 
END 

481 


II     RKGSIPRM1 ,Y, DFRY.NDl*, INLF.FCT.CUTP.AUX) 


1U9 


482 
482 
484 


485 


DIM       .  >  (  1  )  ,DERY  (  1  )  ,AUX(  c,  1)  ,A(4)  ,  B  (  4  )  ,C  (4)  ,PRMT  (  1  ) 
X  =  .'  : f  1(1) 
H*PRKT( i) 


I  0. 
CALL  FCT(X, Y, ULRY) 


487 
486 
489 
4SC 


C 

c 


PREPARATIONS  FOR  RUNCL-KUTTA  METHOO 
2  All).. 


Am  =  .29?R<>  32 
M  \) =1. 707107 
A(4)  =  .16^>6667 


491 
492 


(2>=l. 
(  11*1. 


494 
495 
496 


•  )-2. 
C  (  I  )  =  .  5 
C(2)  =  .2  92  8S)32 


497 
498 


C(  i)=l. 707107 
C(4)=.5 


499 

5CC 


PARATIONS  OF 
DU  5  I=1,NDIM 
AUX( 1,1 )=Y( I ) 


FIRST  RUNGE-KUTTA  STFP 


5C1 
5C2 
5C2 


AUX12, I )=DERY( I ) 
MJX(3,I  )=<  . 
3  AUX(6,1  )=0. 


5C4 


RECORDING  OF  INITIAL  VALUES  OF  THIb  STEP 
7    CALL  OUTP(X,Y,DERY, IREC, NDIM, PRMT )  


5C 


C 

c 


IF( PRMT( 5) )40,8,40 


5C6 
5C7 


STARI  CF  INNERMOST  RUNGE-KUTTA  LOOP 
8  J=l 
10  AJ=A (J) 


(  J) 
CJ=C( J) 

LI  I*1,ND IM 

RL=H*CFRY( I ) 

R2=AJ*(Kl-BJ*AUX(6, I ) ) 

Y(  I  )=Y(  I  )+R2 

♦  ;>  +  R2 
[(6,1) =AUX(6, I )+R2-CJ*Rl 

in?.!1;.!1; 


508 
5C9 
51C 


511 
512 
512 


514 
515 
5U 


11 


If  (  J-4)  1?,  15,  IS 


517 
518 


12  J- J*  I 

IF( J-3)13,14,13 


5  IS        li  X=X  +  .t>:  150 

52C        14  CALL  FCT(X,Y,DERY ) 

521 GOTO  10 

C      END  CF  INNERMOST  RUNGE-KUTTA  LOOP 
C 

C 

522        15  DO  29  I=l,NDIM 

522  AU<( 1, I  )=Y(  I  ) 

524 AUX(2,  I  )=DFRY(1  ) 


525       29  4UXI6, I )=AUX( 3,1 ) 

52fc  CALL  CUT?(X,Y,DERY,  IHLF.NDIM, PRMT) 

527  IF(PQMT(5) KO, 30,40 


5  2 8       30  DO  31  I=ltNDIM 
52S  Y(  I )  =  AUX(  1,  I  ) 

53C       31  DF_^Y(  I  )=AUX(2,  I  ) 


531  GO  TC  8 

532  40  RETURN 

533  ID 
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ABSTRACT 

There  are  many  difficulties  in  using  either  the  classical  multistage 
optimization  techniques  or  dynamic  programming  for  solving  nonlinear 
complex  problems  involving  a  fairly  large  number  of  variables.   The 
former  gives  boundary  value  difficulties  while  the  latter  has  the  dif- 
ficulty of  dimensionality.   The  methods  of  gradients  and  other  techniques 
such  as  quasilinearization  partially  overcome  these  difficulties. 

The  basic  philosophy  of  the  methods  of  gradients  is  fairly  simple. 
First  a  sequence  of  values  of  the  control  vector  is  selected.   Then  the 
gradient  of  the  performance  index  with  respect  to  each  of  the  control 
vector  is  calculated.   Finally  each  control  vector  is  improved  by  moving 
it  in  the  direction  of  the  gradient.   This  improved  sequence  of  control 
vectors  then  becomes  the  basis  for  the  next  iteration. 

The  functional  gradient  technique,  one  of  the  many  versions  of  the 
gradient  methods,  has  been  developed  for  optimal  control  problems.   The 
second  variation  method  overcomes  certain  difficulties  of  the  functional 
gradient  technique.   The  convergence  rate  of  the  second  variation  method, 
provided  the  method  converges,  is  very  fast.   However,  the  initial  guess 
of  the  trajectory  for  the  control  variable  has  to  be  near  the  optimal  tra- 
jectory in  order  to  obtain  convergence.   Too,  the  number  of  equations  to 
be  integrated  and  their  complexity  tend  to  suppress  its  advantage  of  rapid 
convergence. 

First,  the  method  of  second  variation  is  discussed  in  detail.   Then 
the  method  is  applied  to  three  problems  in  the  field  of  production  and 
inventory  control  to  illustrate  the  approach. 


The  first  application  is  a  simple  inventory  model  involving  one  state 
and  one  control  variable.   The  objective  function  is  the  cost 
function,  which  is  to  be  minimized.   The  second  application  is  an  in- 

Qtory  and  advertising  model  where  it  is  desired  to  maximize  the  profit 
function.   This  problem  has  two  state  variables  and  one  control  variable. 
The  last  application  is  that  of  a  chemical  manufacturing  problem  with  ad- 
vertisement.  It  has  six  state  variables  and  three  control  variables. 

These  examples  suggest  that  the  first  variation  method,  of  which  the 
second  variation  method  is  a  natural  evolution,  should  be  used  in  combin- 
ation with  the  second  variation.   The  first  variation  method,  unlike  the 
second  variation,  will  approach  optimum  from  almost  any  realistic  starting 
trajectory.   The  results  of  the  first  variation  method  could  then  be 
used  as  the  starting  trajectories  for  the  second  variation.   In  this  way, 
the  convergence  problem  of  the  second  variation  can  be  partly  overcome. 
Furthermore,  this  combination  provides  a  rapid  convergence  from  almost 
any  realistic  starting  trajectory  for  most  engineering  problems. 


